1 Introduction

In this script, we are building the spatiotemporal prediction model for black carbon in the Denver metro area. This model uses data from 5 sampling campaigns (as detailed in 15_Exploring_BC_Data.R). Two models are developed here: one which uses one temporal basis function (Model A) and one that uses two temporal basis functions (Model B).

Updates to this model are in response to reviewer comments on an earlier manuscript version submitted to ES&T.

A summary table at the end of the document summarizes the performance for each of these models. Overall, Model B performs well and allows us to predict BC concentrations at a weekly temporal resolution from 2009 through 2020.

## Loading required package: Matrix

2 Developing data sets

2.1 Denver data set

Here I’m reading in the filter data sets and getting them into the right shape for the SpatioTemporal package. Here the ST covariates cover the entire study period (2009 - 2020).

data_name <- "Combined_Filter_Data_AEA.csv"
all_data <- read_csv(here::here("Data", data_name))
## Parsed with column specification:
## cols(
##   .default = col_double(),
##   filter_id = col_character(),
##   campaign = col_character(),
##   StartDateTimeLocal = col_date(format = ""),
##   EndDateTimeLocal = col_date(format = ""),
##   sample_week = col_date(format = ""),
##   As_ug_m3 = col_logical(),
##   Cd_ug_m3 = col_logical(),
##   In_ug_m3 = col_logical(),
##   Ni_ug_m3 = col_logical(),
##   Sn_ug_m3 = col_logical(),
##   Te_ug_m3 = col_logical(),
##   st_week = col_date(format = ""),
##   units_pm = col_character(),
##   units_bc = col_character(),
##   units_no2 = col_character(),
##   units_temp = col_character(),
##   units_smoke = col_character(),
##   nn3_bc = col_logical(),
##   area_bc = col_logical(),
##   idw_bc = col_logical()
## )
## See spec(...) for full column specifications.
#' Select a "calibrated" version of the data
#' For now, go with Deming regression-- accounts for variability in the
#' monitor and the UPAS data and the temporal mismatch in TWAs
all_data$bc_ug_m3 <- all_data$bc_ug_m3_dem
all_data$pm_ug_m3 <- all_data$pm_ug_m3_dem

#' List of sites that failed preliminary screening (see 15_Exploring_BC_Data.R)
#' These are all the sites (by campaign) where BC measurements had a CV > 30%
#' Note that all of these are in Campaign 4
cv_drop <- read_csv(here::here("Data/Dropped_Sites_by_Campaign.csv"))
## Parsed with column specification:
## cols(
##   site_id = col_double(),
##   campaign = col_character(),
##   n_filters = col_double(),
##   mean = col_double(),
##   sd = col_double(),
##   cv = col_double(),
##   median = col_double(),
##   IQR = col_double(),
##   range = col_double(),
##   cv_gt_50 = col_double(),
##   cv_gt_30 = col_double()
## )
filter_data <- filter(all_data, is.na(filter_id) | filter_id != "080310027") %>%
  filter(is.na(indoor) | indoor == 0) %>%
  filter(is.na(is_blank) | is_blank == 0) %>% 
  #' QA filters
  filter(is.na(bc_below_lod) | bc_below_lod == 0) %>% 
  filter(is.na(negative_pm_mass) | negative_pm_mass == 0) %>% 
  filter(is.na(potential_contamination) | potential_contamination == 0)

dist_data <- bind_rows(filter(filter_data, is.na(campaign) | campaign != "Campaign4"),
                       filter(filter_data, campaign == "Campaign4" & !(site_id %in% cv_drop$site_id))) %>%
  arrange(st_week)
head(dist_data$sample_week)
## [1] NA NA NA NA NA NA
tail(dist_data$sample_week)
## [1] NA NA NA NA NA NA
head(dist_data$st_week)
## [1] "2008-12-29" "2008-12-29" "2008-12-29" "2008-12-29" "2008-12-29"
## [6] "2008-12-29"
tail(dist_data$st_week)
## [1] "2020-12-28" "2020-12-28" "2020-12-28" "2020-12-28" "2020-12-28"
## [6] "2020-12-28"
length(unique(dist_data$filter_id))
## [1] 848
summary(dist_data$bc_ug_m3)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
##    0.51    0.99    1.17    1.27    1.46    2.27   41600
quantile(dist_data$bc_ug_m3, probs = c(0.05, 0.95), na.rm = T)
##        5%       95% 
## 0.8689785 2.0648167
sd(dist_data$bc_ug_m3, na.rm = T)
## [1] 0.3498975
sd(dist_data$bc_ug_m3, na.rm = T) / mean(dist_data$bc_ug_m3, na.rm = T)
## [1] 0.2764677
#' All of the dropped data should be in Campaign 4
dropped_data <- filter(filter_data, campaign == "Campaign4" & site_id %in% cv_drop$site_id)
length(unique(dropped_data$filter_id))
## [1] 97
summary(dropped_data$bc_ug_m3)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
## -2.7752  0.6319  1.4887  1.1935  1.9480  2.1754
quantile(dropped_data$bc_ug_m3, probs = c(0.05, 0.95))
##         5%        95% 
## -0.3137664  2.1225874
sd(dropped_data$bc_ug_m3)
## [1] 0.9202973
sd(dropped_data$bc_ug_m3) / mean(dropped_data$bc_ug_m3)
## [1] 0.7711097
#' central site data
central_data <- filter(all_data, filter_id == "080310027") %>%
  arrange(st_week)
head(central_data$sample_week)
## [1] "2008-12-29" "2009-01-05" "2009-01-12" "2009-01-19" "2009-01-26"
## [6] "2009-02-02"
tail(central_data$sample_week)
## [1] "2020-11-23" "2020-11-30" "2020-12-07" "2020-12-14" "2020-12-21"
## [6] "2020-12-28"
head(central_data$st_week)
## [1] "2008-12-29" "2009-01-05" "2009-01-12" "2009-01-19" "2009-01-26"
## [6] "2009-02-02"
tail(central_data$st_week)
## [1] "2020-11-23" "2020-11-30" "2020-12-07" "2020-12-14" "2020-12-21"
## [6] "2020-12-28"
#' Add prefix to site_ids
unique(dist_data$site_id)
##  [1]  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25
## [26] 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50
## [51] 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68
dist_data$site_id <- paste0("d_", dist_data$site_id)
unique(dist_data$site_id)
##  [1] "d_1"  "d_2"  "d_3"  "d_4"  "d_5"  "d_6"  "d_7"  "d_8"  "d_9"  "d_10"
## [11] "d_11" "d_12" "d_13" "d_14" "d_15" "d_16" "d_17" "d_18" "d_19" "d_20"
## [21] "d_21" "d_22" "d_23" "d_24" "d_25" "d_26" "d_27" "d_28" "d_29" "d_30"
## [31] "d_31" "d_32" "d_33" "d_34" "d_35" "d_36" "d_37" "d_38" "d_39" "d_40"
## [41] "d_41" "d_42" "d_43" "d_44" "d_45" "d_46" "d_47" "d_48" "d_49" "d_50"
## [51] "d_51" "d_52" "d_53" "d_54" "d_55" "d_56" "d_57" "d_58" "d_59" "d_60"
## [61] "d_61" "d_62" "d_63" "d_64" "d_65" "d_66" "d_67" "d_68"
unique(central_data$site_id)
## [1] 16
central_data$site_id <- "central"
unique(central_data$site_id)
## [1] "central"
#' put these data sets back together
lur_data <- bind_rows(dist_data, central_data) %>%
  arrange(site_id, st_week)
head(lur_data$st_week)
## [1] "2008-12-29" "2009-01-05" "2009-01-12" "2009-01-19" "2009-01-26"
## [6] "2009-02-02"
tail(lur_data$st_week)
## [1] "2020-11-23" "2020-11-30" "2020-12-07" "2020-12-14" "2020-12-21"
## [6] "2020-12-28"
#' Log-transformed BC observations with NA values for the missing dates
#' Note that campaign 5 has duplicate measures, so we need to average them
bc_obs <- lur_data %>%
  dplyr::select(site_id, st_week, bc_ug_m3) %>%
  group_by(site_id, st_week) %>%
  summarize(bc_ug_m3 = mean(bc_ug_m3, na.rm = T)) %>%
  mutate(bc_ug_m3 = ifelse(is.nan(bc_ug_m3), NA, bc_ug_m3)) %>%
  mutate(log_bc = log(bc_ug_m3)) %>%
  dplyr::select(-bc_ug_m3) 

#' Pivot to a wide data frame
bc_obs2 <- bc_obs %>%
  pivot_wider(id_cols = st_week, names_from = site_id, values_from = log_bc) %>%
  # names_from = site_id, values_from = bc_ug_m3) %>%
  as.data.frame() %>%
  arrange(st_week)
rownames(bc_obs2) <- bc_obs2$st_week
bc_obs2$st_week <- NULL
bc_obs2 <- as.matrix(bc_obs2)

bc_weeks <- rownames(bc_obs2)
tail(bc_weeks)
## [1] "2020-11-23" "2020-11-30" "2020-12-07" "2020-12-14" "2020-12-21"
## [6] "2020-12-28"
#' Spatial covariates (scaled)
load(here::here("Results", "BC_LASSO_Model_5C.Rdata"))

lasso_coef_df <- data.frame(name = log_bc_lasso_coef4@Dimnames[[1]][log_bc_lasso_coef4@i + 1],
                            coefficient = log_bc_lasso_coef4@x)
covars <- as.character(lasso_coef_df$name[-1])
covars 
##  [1] "tree_cover_100"  "impervious_2500" "low_int_100"     "med_int_50"     
##  [5] "med_int_2500"    "high_int_100"    "pop_den_50"      "dist_m_cafos"   
##  [9] "dist_m_og_wells" "dist_m_wwtp"     "aadt_100"        "aadt_250"       
## [13] "aadt_1000"
covar_fun <- paste("~", paste(covars, collapse = " + "))
covar_fun
## [1] "~ tree_cover_100 + impervious_2500 + low_int_100 + med_int_50 + med_int_2500 + high_int_100 + pop_den_50 + dist_m_cafos + dist_m_og_wells + dist_m_wwtp + aadt_100 + aadt_250 + aadt_1000"
bc_sp_cov <- dplyr::select(lur_data, site_id, lon, lat, all_of(covars)) %>%
  distinct() %>%
  mutate_at(.vars = vars(all_of(covars)),
            scale)

bc_sp_cov <- bc_sp_cov %>%
  rename(ID = site_id) %>%
  mutate(lon2 = lon, lat2 = lat) %>%
  st_as_sf(coords = c("lon2", "lat2"), crs = ll_wgs84) %>%
  st_transform(crs = albers)
sp_coords <- do.call(rbind, st_geometry(bc_sp_cov)) %>% as_tibble()
names(sp_coords) <- c("x", "y")

bc_sp_cov <- bind_cols(bc_sp_cov, sp_coords) %>%
  st_set_geometry(NULL) %>%
  as.data.frame()

bc_sp_cov$type <- ifelse(bc_sp_cov$ID == "central", "central", "dist")
bc_sp_cov$type <- as.factor(bc_sp_cov$type)

cor(bc_sp_cov[,covars])
##                 tree_cover_100 impervious_2500  low_int_100  med_int_50
## tree_cover_100     1.000000000      0.05592620  0.227568688 -0.41656599
## impervious_2500    0.055926204      1.00000000  0.046627086 -0.13878441
## low_int_100        0.227568688      0.04662709  1.000000000 -0.11614316
## med_int_50        -0.416565993     -0.13878441 -0.116143165  1.00000000
## med_int_2500       0.005240397      0.89874242  0.172586606 -0.13029847
## high_int_100      -0.335160074      0.48006634 -0.444088483  0.11084350
## pop_den_50         0.201795804      0.39972348  0.483551971  0.05903168
## dist_m_cafos       0.262782544     -0.02217726 -0.161007912 -0.09924381
## dist_m_og_wells    0.304301808     -0.11331048  0.008034836 -0.21284246
## dist_m_wwtp        0.140092353     -0.49234185 -0.075238446  0.02725541
## aadt_100          -0.176079221      0.40911376 -0.255918459  0.06413761
## aadt_250          -0.007357196      0.42742607 -0.261080670  0.03197137
## aadt_1000          0.085749987      0.43050282 -0.112576947 -0.19807436
##                 med_int_2500 high_int_100  pop_den_50 dist_m_cafos
## tree_cover_100   0.005240397  -0.33516007  0.20179580   0.26278254
## impervious_2500  0.898742424   0.48006634  0.39972348  -0.02217726
## low_int_100      0.172586606  -0.44408848  0.48355197  -0.16100791
## med_int_50      -0.130298466   0.11084350  0.05903168  -0.09924381
## med_int_2500     1.000000000   0.39325577  0.51099015  -0.24164427
## high_int_100     0.393255768   1.00000000 -0.10833432  -0.01905821
## pop_den_50       0.510990152  -0.10833432  1.00000000  -0.11900734
## dist_m_cafos    -0.241644265  -0.01905821 -0.11900734   1.00000000
## dist_m_og_wells -0.251552596  -0.21708310 -0.02501575   0.61986048
## dist_m_wwtp     -0.456542801  -0.30637142 -0.14260842   0.15352999
## aadt_100         0.261463299   0.71529092 -0.08808862   0.10155682
## aadt_250         0.276833428   0.60158329 -0.10795874   0.23112012
## aadt_1000        0.349181488   0.40408074  0.01155408   0.12765618
##                 dist_m_og_wells dist_m_wwtp    aadt_100     aadt_250
## tree_cover_100      0.304301808  0.14009235 -0.17607922 -0.007357196
## impervious_2500    -0.113310476 -0.49234185  0.40911376  0.427426070
## low_int_100         0.008034836 -0.07523845 -0.25591846 -0.261080670
## med_int_50         -0.212842464  0.02725541  0.06413761  0.031971369
## med_int_2500       -0.251552596 -0.45654280  0.26146330  0.276833428
## high_int_100       -0.217083104 -0.30637142  0.71529092  0.601583294
## pop_den_50         -0.025015754 -0.14260842 -0.08808862 -0.107958743
## dist_m_cafos        0.619860480  0.15352999  0.10155682  0.231120116
## dist_m_og_wells     1.000000000  0.35073325 -0.06090070  0.109442838
## dist_m_wwtp         0.350733252  1.00000000 -0.20875176 -0.024182445
## aadt_100           -0.060900696 -0.20875176  1.00000000  0.806259695
## aadt_250            0.109442838 -0.02418244  0.80625970  1.000000000
## aadt_1000           0.091175296 -0.05838643  0.50962569  0.635217540
##                   aadt_1000
## tree_cover_100   0.08574999
## impervious_2500  0.43050282
## low_int_100     -0.11257695
## med_int_50      -0.19807436
## med_int_2500     0.34918149
## high_int_100     0.40408074
## pop_den_50       0.01155408
## dist_m_cafos     0.12765618
## dist_m_og_wells  0.09117530
## dist_m_wwtp     -0.05838643
## aadt_100         0.50962569
## aadt_250         0.63521754
## aadt_1000        1.00000000
#' NO2, and smoke spatiotemporal predictors
#' NO2 at each site estimated using IDW
bc_st_no2 <- dplyr::select(lur_data, site_id, st_week, idw_no2) %>%
  distinct() %>%
  arrange(st_week) %>%
  pivot_wider(id_cols = st_week,
              names_from = site_id, values_from = idw_no2) %>%
  as.data.frame()
rownames(bc_st_no2) <- bc_st_no2$st_week
bc_st_no2$st_week <- NULL
bc_st_no2 <- as.matrix(bc_st_no2)

no2_weeks <- rownames(bc_st_no2)
tail(no2_weeks)
## [1] "2020-11-23" "2020-11-30" "2020-12-07" "2020-12-14" "2020-12-21"
## [6] "2020-12-28"
setdiff(bc_weeks, no2_weeks)
## character(0)
#' Smoke day indicator variable based on WFS paper method (see Brey and Fischer, 2016)
#' based on any smoke variable in the area using the 2 sd cut-off (area_smoke_2sd)
bc_st_smk <- dplyr::select(lur_data, site_id, st_week, area_smoke_2sd) %>%
  distinct() %>%
  arrange(st_week) %>%
  pivot_wider(id_cols = st_week,
              names_from = site_id, values_from = area_smoke_2sd) %>%
  as.data.frame()
rownames(bc_st_smk) <- bc_st_smk$st_week
bc_st_smk$st_week <- NULL
bc_st_smk <- as.matrix(bc_st_smk)

smk_weeks <- rownames(bc_st_smk)
tail(smk_weeks)
## [1] "2020-11-23" "2020-11-30" "2020-12-07" "2020-12-14" "2020-12-21"
## [6] "2020-12-28"
setdiff(bc_weeks, smk_weeks)
## character(0)

Create a new denver.data object.

denver.data <- createSTdata(obs = bc_obs2,
                            covars = bc_sp_cov,
                            SpatioTemporal = list(bc_st_no2 = bc_st_no2,
                                                  bc_st_smk = bc_st_smk))
D <- createDataMatrix(denver.data)
denver.data
## STdata-object with:
##  No. locations: 69 (observed: 64)
##  No. time points: 231 (observed: 231)
##  No. obs: 1021
## 
## Only constant temporal trend,with dates:
##  2016-02-08 to 2020-11-30
## 
## 19 covariate(s):
##  [1] "ID"              "lon"             "lat"             "tree_cover_100" 
##  [5] "impervious_2500" "low_int_100"     "med_int_50"      "med_int_2500"   
##  [9] "high_int_100"    "pop_den_50"      "dist_m_cafos"    "dist_m_og_wells"
## [13] "dist_m_wwtp"     "aadt_100"        "aadt_250"        "aadt_1000"      
## [17] "x"               "y"               "type"           
## 
## 2 spatio-temporal covariate(s):
## [1] "bc_st_no2" "bc_st_smk"
## 
## All sites:
## central    dist 
##       1      68 
## Observed:
## central    dist 
##       1      63 
## 
## For central:
##   Number of obs: 231
##   Dates: 2016-02-08 to 2020-11-30
## For dist:
##   Number of obs: 790
##   Dates: 2018-05-07 to 2020-04-06

2.2 Creating the NO2, BC, and PM data objects

Next we need to combine monitoring data for NO2, BC, and PM2.5 in order to calculate the time trends.

First up: NO2

Plot the time trends for NO2 at the monitors

Original units (ppb)

no2_ts <- ggplot(data = no2_data, aes(x = Date_Local, y = Arithmetic_Mean)) +
  geom_point(aes(color = as.factor(monitor_id)), shape = 20, size = 1) +
  geom_smooth(se = F, color = "black", method = "gam") +
  scale_color_viridis(name = "Monitor ID", discrete = T) +
  scale_x_date(date_breaks = "6 months", date_labels = "%m-%Y") +
  facet_wrap(. ~ monitor_id, scales = "fixed", ncol = 2) +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) +
  simple_theme
no2_ts
## `geom_smooth()` using formula 'y ~ s(x, bs = "cs")'

Log-transformed NO2

no2_ts2 <- ggplot(data = no2_data, aes(x = Date_Local, y = log(Arithmetic_Mean))) +
  geom_point(aes(color = as.factor(monitor_id)), shape = 20, size = 1) +
  geom_smooth(se = F, color = "black", method = "gam") +
  scale_color_viridis(name = "Monitor ID", discrete = T) +
  scale_x_date(date_breaks = "6 months", date_labels = "%m-%Y") +
  facet_wrap(. ~ monitor_id, scales = "fixed", ncol = 2) +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) +
  simple_theme
no2_ts2
## `geom_smooth()` using formula 'y ~ s(x, bs = "cs")'

Organize the observations and spatial data. Use the log-transformed NO2 here

#' NO2 observations
no2_obs <- no2_data %>%
  st_set_geometry(NULL) %>%
  dplyr::select(monitor_id, Date_Local, Arithmetic_Mean) %>%
  mutate(st_week = as.character(as.Date(cut(as.Date(Date_Local), "week")))) %>%
  filter(st_week %in% bc_weeks) %>%
  group_by(monitor_id, st_week) %>%
  mutate(log_mean = log(Arithmetic_Mean)) %>%
  summarize(no2 = mean(log_mean, na.rm=T)) %>%
  pivot_wider(id_cols = st_week,
              names_from = monitor_id, values_from = no2) %>%
  as.data.frame() %>%
  arrange(st_week)
head(no2_obs)
head(no2_obs$st_week)
## [1] "2008-12-29" "2009-01-05" "2009-01-12" "2009-01-19" "2009-01-26"
## [6] "2009-02-02"
tail(no2_obs$st_week)
## [1] "2020-11-02" "2020-11-09" "2020-11-16" "2020-11-23" "2020-11-30"
## [6] "2020-12-07"
# NO2 SP object
no2_data2 <- read_csv(here::here("Data", "Monitor_NO2_Data_AEA.csv")) %>%
  mutate(year = year(Date_Local)) %>%
  dplyr::select(-year) %>%
  st_as_sf(wkt = "WKT", crs = albers) %>%
  mutate(monitor_id = paste0(monitor_id, "_no2")) %>%
  filter(monitor_id %in% colnames(no2_obs))
## Parsed with column specification:
## cols(
##   .default = col_character(),
##   Parameter_Code = col_double(),
##   POC = col_double(),
##   Date_Local = col_date(format = ""),
##   Observation_Count = col_double(),
##   Observation_Percent = col_double(),
##   Arithmetic_Mean = col_double(),
##   Max_Value = col_double(),
##   `1st_Max_Hour` = col_double(),
##   AQI = col_double(),
##   Date_of_Last_Change = col_date(format = "")
## )
## See spec(...) for full column specifications.
no2_coords <- do.call(rbind, st_geometry(no2_data2)) %>% as_tibble()
names(no2_coords) <- c("x", "y")
no2_sp_lonlat <- no2_data2 %>%
  st_transform(crs = ll_wgs84)
no2_coords2 <- do.call(rbind, st_geometry(no2_sp_lonlat)) %>%
  as_tibble() %>% setNames(c("lon","lat"))

no2_sp_cov <- st_set_geometry(no2_data2, NULL) %>%
  dplyr::select(monitor_id) %>%
  rename(ID = monitor_id)
no2_sp_cov <- bind_cols(no2_sp_cov, no2_coords, no2_coords2) %>%
  as.data.frame() %>%
  distinct()

Next, BC at the central site monitor

#' BC central site concentrations
bc_cent_data <- read_csv(here::here("Data", "Monitor_BC_Data_AEA.csv")) %>%
  filter(!is.na(Arithmetic_Mean)) %>%
  st_as_sf(wkt = "WKT", crs = albers) %>%
  mutate(County_Code = str_sub(monitor_id, start = 3, end = 5)) %>%
  filter(County_Code %in% counties) %>%
  arrange(Date_Local, monitor_id) %>%
  mutate(Arithmetic_Mean = ifelse(Arithmetic_Mean <= 0.005, NA, Arithmetic_Mean))
## Parsed with column specification:
## cols(
##   .default = col_character(),
##   Parameter_Code = col_double(),
##   POC = col_double(),
##   Pollutant_Standard = col_logical(),
##   Date_Local = col_date(format = ""),
##   Observation_Count = col_double(),
##   Observation_Percent = col_double(),
##   Arithmetic_Mean = col_double(),
##   Max_Value = col_double(),
##   `1st_Max_Hour` = col_double(),
##   AQI = col_logical(),
##   Method_Code = col_double(),
##   Date_of_Last_Change = col_date(format = "")
## )
## See spec(...) for full column specifications.
head(bc_cent_data$Date_Local)
## [1] "2016-02-08" "2016-02-09" "2016-02-10" "2016-02-11" "2016-02-12"
## [6] "2016-02-13"
tail(bc_cent_data$Date_Local)
## [1] "2020-11-25" "2020-11-26" "2020-11-27" "2020-11-28" "2020-11-29"
## [6] "2020-11-30"
#' Add an identifier to differentiate these measurements from collocated pollutants at the same site
bc_cent_data$monitor_id <- paste0(bc_cent_data$monitor_id, "_bc")
unique(bc_cent_data$monitor_id)
## [1] "080310027_bc"

Plot the time trends for BC at the central site monitor

Original units (ug/m3)

bc_ts <- ggplot(data = bc_cent_data, aes(x = Date_Local, y = Arithmetic_Mean)) +
  geom_point(aes(color = as.factor(monitor_id)), shape = 20, size = 1) +
  geom_smooth(se = F, color = "black", method = "gam") +
  scale_color_viridis(name = "Monitor ID", discrete = T) +
  scale_x_date(date_breaks = "6 months", date_labels = "%m-%Y") +
  #facet_wrap(. ~ monitor_id, scales = "fixed", ncol = 2) +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) +
  simple_theme
bc_ts
## `geom_smooth()` using formula 'y ~ s(x, bs = "cs")'

Log-transformed BC

bc_ts2 <- ggplot(data = bc_cent_data, aes(x = Date_Local, y = log(Arithmetic_Mean))) +
  geom_point(aes(color = as.factor(monitor_id)), shape = 20, size = 1) +
  geom_smooth(se = F, color = "black") +
  scale_color_viridis(name = "Monitor ID", discrete = T) +
  scale_x_date(date_breaks = "6 months", date_labels = "%m-%Y") +
  #facet_wrap(. ~ monitor_id, scales = "fixed", ncol = 2) +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) +
  simple_theme
bc_ts2
## `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'

Organize the observations and spatial data

Use log-transformed data here

#' BC observations
bc_cent_obs <- bc_cent_data %>%
  st_set_geometry(NULL) %>%
  dplyr::select(monitor_id, Date_Local, Arithmetic_Mean) %>%
  mutate(st_week = as.character(as.Date(cut(as.Date(Date_Local), "week")))) %>%
  filter(st_week %in% bc_weeks) %>%
  group_by(monitor_id, st_week) %>%
  mutate(log_mean = log(Arithmetic_Mean)) %>%
  summarize(bc = mean(log_mean, na.rm=T)) %>%
  pivot_wider(id_cols = st_week,
              names_from = monitor_id, values_from = bc) %>%
  as.data.frame() %>%
  arrange(st_week)
head(bc_cent_obs)
head(bc_cent_obs$st_week)
## [1] "2016-02-08" "2016-02-15" "2016-02-22" "2016-02-29" "2016-03-07"
## [6] "2016-03-14"
tail(bc_cent_obs$st_week)
## [1] "2020-10-26" "2020-11-02" "2020-11-09" "2020-11-16" "2020-11-23"
## [6] "2020-11-30"
# BC Central Site SP object
bc_cent_data2 <- read_csv(here::here("Data", "Monitor_BC_Data_AEA.csv")) %>%
  mutate(year = year(Date_Local)) %>%
  dplyr::select(-year) %>%
  st_as_sf(wkt = "WKT", crs = albers) %>%
  mutate(monitor_id = paste0(monitor_id, "_bc")) %>%
  filter(monitor_id %in% colnames(bc_cent_obs))
## Parsed with column specification:
## cols(
##   .default = col_character(),
##   Parameter_Code = col_double(),
##   POC = col_double(),
##   Pollutant_Standard = col_logical(),
##   Date_Local = col_date(format = ""),
##   Observation_Count = col_double(),
##   Observation_Percent = col_double(),
##   Arithmetic_Mean = col_double(),
##   Max_Value = col_double(),
##   `1st_Max_Hour` = col_double(),
##   AQI = col_logical(),
##   Method_Code = col_double(),
##   Date_of_Last_Change = col_date(format = "")
## )
## See spec(...) for full column specifications.
bc_cent_coords <- do.call(rbind, st_geometry(bc_cent_data2)) %>% as_tibble()
names(bc_cent_coords) <- c("x", "y")
bc_cent_sp_lonlat <- bc_cent_data2 %>%
  st_transform(crs = ll_wgs84)
bc_cent_coords2 <- do.call(rbind, st_geometry(bc_cent_sp_lonlat)) %>%
  as_tibble() %>% setNames(c("lon","lat"))

bc_cent_sp_cov <- st_set_geometry(bc_cent_data2, NULL) %>%
  dplyr::select(monitor_id) %>%
  rename(ID = monitor_id)
bc_cent_sp_cov <- bind_cols(bc_cent_sp_cov, bc_cent_coords, bc_cent_coords2) %>%
  as.data.frame() %>%
  distinct()

Lastly, PM2.5 from select monitors with long-term records

Plot the time trends for PM2.5 at the monitors

Original units (ug/m3)

pm_ts <- ggplot(data = pm_data, aes(x = Date_Local, y = Arithmetic_Mean)) +
  geom_point(aes(color = as.factor(monitor_id)), shape = 20, size = 1) +
  geom_smooth(se = F, color = "black") +
  scale_color_viridis(name = "Monitor ID", discrete = T) +
  scale_x_date(date_breaks = "6 months", date_labels = "%m-%Y") +
  facet_wrap(. ~ monitor_id, scales = "fixed", ncol = 2) +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) +
  simple_theme
pm_ts
## `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'

Log-transformed PM2.5

pm_ts2 <- ggplot(data = pm_data, aes(x = Date_Local, y = log(Arithmetic_Mean))) +
  geom_point(aes(color = as.factor(monitor_id)), shape = 20, size = 1) +
  geom_smooth(se = F, color = "black", method = "gam") +
  scale_color_viridis(name = "Monitor ID", discrete = T) +
  scale_x_date(date_breaks = "6 months", date_labels = "%m-%Y") +
  facet_wrap(. ~ monitor_id, scales = "fixed", ncol = 2) +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) +
  simple_theme
pm_ts2
## `geom_smooth()` using formula 'y ~ s(x, bs = "cs")'

Organize the observations and spatial data

Use the log-transformed data here

#' PM2.5 observations
pm_obs <- pm_data %>%
  st_set_geometry(NULL) %>%
  dplyr::select(monitor_id, Date_Local, Arithmetic_Mean) %>%
  mutate(st_week = as.character(as.Date(cut(as.Date(Date_Local), "week")))) %>%
  filter(st_week %in% bc_weeks) %>%
  group_by(monitor_id, st_week) %>%
  mutate(log_mean = log(Arithmetic_Mean)) %>%
  summarize(pm = mean(log_mean, na.rm=T)) %>%
  pivot_wider(id_cols = st_week,
              names_from = monitor_id, values_from = pm) %>%
  as.data.frame() %>%
  arrange(st_week)
head(pm_obs)
head(pm_obs$st_week)
## [1] "2008-12-29" "2009-01-05" "2009-01-12" "2009-01-19" "2009-01-26"
## [6] "2009-02-02"
tail(pm_obs$st_week)
## [1] "2020-10-26" "2020-11-02" "2020-11-09" "2020-11-16" "2020-11-23"
## [6] "2020-11-30"
# PM SP object
pm_data2 <- read_csv(here::here("Data", "Monitor_PM_Data_AEA.csv")) %>%
  mutate(year = year(Date_Local)) %>%
  dplyr::select(-year) %>%
  st_as_sf(wkt = "WKT", crs = albers) %>%
  mutate(monitor_id = paste0(monitor_id, "_pm")) %>%
  filter(monitor_id %in% colnames(pm_obs))
## Parsed with column specification:
## cols(
##   .default = col_character(),
##   Parameter_Code = col_double(),
##   POC = col_double(),
##   Date_Local = col_date(format = ""),
##   Observation_Count = col_double(),
##   Observation_Percent = col_double(),
##   Arithmetic_Mean = col_double(),
##   Max_Value = col_double(),
##   `1st_Max_Hour` = col_double(),
##   AQI = col_double(),
##   Method_Code = col_double(),
##   Date_of_Last_Change = col_date(format = "")
## )
## See spec(...) for full column specifications.
pm_coords <- do.call(rbind, st_geometry(pm_data2)) %>% as_tibble()
names(pm_coords) <- c("x", "y")
pm_sp_lonlat <- pm_data2 %>%
  st_transform(crs = ll_wgs84)
pm_coords2 <- do.call(rbind, st_geometry(pm_sp_lonlat)) %>%
  as_tibble() %>% setNames(c("lon","lat"))

pm_sp_cov <- st_set_geometry(pm_data2, NULL) %>%
  dplyr::select(monitor_id) %>%
  rename(ID = monitor_id)
pm_sp_cov <- bind_cols(pm_sp_cov, pm_coords, pm_coords2) %>%
  as.data.frame() %>%
  distinct()

Combine the observations and the spatial covariates Scale the pollutant measurements

#' Join log-transformed observations
pol_obs <- left_join(no2_obs, bc_cent_obs, by = "st_week") %>%
  left_join(pm_obs, by = "st_week")
glimpse(pol_obs)
## Rows: 622
## Columns: 18
## $ st_week         <chr> "2008-12-29", "2009-01-05", "2009-01-12", "2009-01-19…
## $ `080013001_no2` <dbl> 2.4313422, 2.3337153, 3.0953820, 2.4960668, 2.6445846…
## $ `080310002_no2` <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, N…
## $ `080310026_no2` <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, N…
## $ `080310027_no2` <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, N…
## $ `080310028_no2` <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, N…
## $ `080310027_bc`  <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, N…
## $ `080010006_pm`  <dbl> 1.394118, 1.247789, 1.934844, 2.298456, 1.598321, 1.7…
## $ `080010008_pm`  <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, N…
## $ `080050005_pm`  <dbl> 1.2490759, 0.6817687, 1.3697744, 1.7730816, 1.1939225…
## $ `080130003_pm`  <dbl> 0.9679299, 0.8243293, 2.1047287, 1.7212142, 1.9059908…
## $ `080130012_pm`  <dbl> 1.0919008, 0.8908546, 1.2511276, 1.2703657, 0.8047190…
## $ `080310002_pm`  <dbl> 1.335209, 1.188381, 1.767997, 2.027844, 1.486213, 1.7…
## $ `080310023_pm`  <dbl> 1.661197, 1.513911, 1.919262, 1.996528, 1.909773, 1.7…
## $ `080310025_pm`  <dbl> 1.1537863, 0.9359011, 1.3262687, 1.9969858, 1.6521301…
## $ `080310026_pm`  <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, N…
## $ `080310027_pm`  <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, N…
## $ `080310028_pm`  <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, N…
rownames(pol_obs) <- pol_obs$st_week
pol_obs$st_week <- NULL
pol_obs <- as.matrix(pol_obs)

#' Scale the measurements to avoid issues where units are different
#' Remember! These have been log-transformed before scaling
#' Scale the data by columns (monitors)
rnames <- rownames(pol_obs)
cnames <- colnames(pol_obs)

pol_obs <- apply(pol_obs, 2, scale)

rownames(pol_obs) <- rnames
colnames(pol_obs) <- cnames

head(rownames(pol_obs))
## [1] "2008-12-29" "2009-01-05" "2009-01-12" "2009-01-19" "2009-01-26"
## [6] "2009-02-02"
tail(rownames(pol_obs))
## [1] "2020-11-02" "2020-11-09" "2020-11-16" "2020-11-23" "2020-11-30"
## [6] "2020-12-07"
head(colnames(pol_obs))
## [1] "080013001_no2" "080310002_no2" "080310026_no2" "080310027_no2"
## [5] "080310028_no2" "080310027_bc"
class(pol_obs)
## [1] "matrix" "array"
dim(pol_obs)
## [1] 622  17
#' SP variables
pol_sp_cov <- bind_rows(no2_sp_cov, bc_cent_sp_cov, pm_sp_cov)
pol_sp_cov

Create the data object

pol_STdata <- createSTdata(obs = pol_obs,
                           covars = pol_sp_cov)
print(pol_STdata)
## STdata-object with:
##  No. locations: 17 (observed: 16)
##  No. time points: 622 (observed: 622)
##  No. obs: 6560
## 
## Only constant temporal trend,with dates:
##  2008-12-29 to 2020-12-07
## 
## 5 covariate(s):
## [1] "ID"  "x"   "y"   "lon" "lat"
## 
## No spatio-temporal covariates.

2.3 Generate the temporal trend

The cross validation results for generating the time trends suggest 2 basis functions is the way to go, but let’s see if that holds up with our data

D_pol <- createDataMatrix(pol_STdata)
#D_pol

n_years <- length(2009:2020)
n_years*4
## [1] 48
n_years*8
## [1] 96
#' Trying 4 df per year
pol.SVD.cv.4py <- SVDsmoothCV(D_pol, 0:4, df = n_years * 4)
print(pol.SVD.cv.4py)
## Result of SVDsmoothCV, average of CV-statistics:
##                 MSE        R2          AIC         BIC
## n.basis.0 0.9975610 0.0000000    0.9987781    5.095695
## n.basis.1 0.8013682 0.1967023  -96.6124965  -88.418663
## n.basis.2 0.6358644 0.3625750 -225.3897791 -213.099029
## n.basis.3 0.6297167 0.3687387 -229.4297802 -213.042114
## n.basis.4 0.6084604 0.3900698 -240.5027300 -220.018147
plot(pol.SVD.cv.4py)

#' Trying 8 df per year
pol.SVD.cv.8py <- SVDsmoothCV(D_pol, 0:4, df = n_years * 8)
print(pol.SVD.cv.8py)
## Result of SVDsmoothCV, average of CV-statistics:
##                 MSE        R2          AIC         BIC
## n.basis.0 0.9975610 0.0000000    0.9987781    5.095695
## n.basis.1 0.7630671 0.2350911 -120.2211615 -112.027328
## n.basis.2 0.5825979 0.4159654 -265.3105795 -253.019830
## n.basis.3 0.5739736 0.4246129 -271.3027690 -254.915103
## n.basis.4 0.5615078 0.4371173 -278.8000243 -258.315442
plot(pol.SVD.cv.8py)

2.5 Update the Denver data set

The next step is to update the trend data for the denver.data object and see if these trend data fit our observations.

Comparing the plots using one basis function vs. two basis functions, it looks like using two basis functions might be better than 1. There was some residual autocorrelation in the central site monitor data when using just one basis function (see the fourth plot in the first series of plots below). Using 8 degrees of freedom did not help.

Using two basis functions helped somewhat to reduce this autocorrelation, but it still persists. There are still some noticable trends in the residuals for the central monitoring site

2.5.1 1 basis function, 4 df per year

First update the trend data, then plot the trends for three distributed sites and the central BC monitoring site.

denver.data2.1 <- denver.data
denver.data2.1$trend <- pol_STdata1$trend
print(denver.data2.1)
## STdata-object with:
##  No. locations: 69 (observed: 64)
##  No. time points: 622 (observed: 231)
##  No. obs: 1021
## 
## Trend with 1 basis function(s):
## [1] "V1"
## with dates:
##  2008-12-29 to 2020-12-07
## 
## 19 covariate(s):
##  [1] "ID"              "lon"             "lat"             "tree_cover_100" 
##  [5] "impervious_2500" "low_int_100"     "med_int_50"      "med_int_2500"   
##  [9] "high_int_100"    "pop_den_50"      "dist_m_cafos"    "dist_m_og_wells"
## [13] "dist_m_wwtp"     "aadt_100"        "aadt_250"        "aadt_1000"      
## [17] "x"               "y"               "type"           
## 
## 2 spatio-temporal covariate(s):
## [1] "bc_st_no2" "bc_st_smk"
## 
## All sites:
## central    dist 
##       1      68 
## Observed:
## central    dist 
##       1      63 
## 
## For central:
##   Number of obs: 231
##   Dates: 2016-02-08 to 2020-11-30
## For dist:
##   Number of obs: 790
##   Dates: 2018-05-07 to 2020-04-06
par(mar=c(3.3,3.3,1.5,1), mgp=c(2,1,0))
layout(matrix(c(1,1,2,2,3,4), 3, 2, byrow=TRUE))
plot(denver.data2.1, "obs", ID="d_1", xlab="", ylab="BC (log ug/m3)",
     main="Temporal trend d_1")
plot(denver.data2.1, "res", ID="d_1", xlab="", ylab="BC (log ug/m3)")
plot(denver.data2.1, "acf", ID="d_1")
plot(denver.data2.1, "pacf", ID="d_1")

par(mar=c(3.3,3.3,1.5,1), mgp=c(2,1,0))
layout(matrix(c(1,1,2,2,3,4), 3, 2, byrow=TRUE))
plot(denver.data2.1, "obs", ID="d_10", xlab="", ylab="BC (log ug/m3)",
     main="Temporal trend d_43")
plot(denver.data2.1, "res", ID="d_10", xlab="", ylab="BC (log ug/m3)")
plot(denver.data2.1, "acf", ID="d_10")
plot(denver.data2.1, "pacf", ID="d_10")

par(mar=c(3.3,3.3,1.5,1), mgp=c(2,1,0))
layout(matrix(c(1,1,2,2,3,4), 3, 2, byrow=TRUE))
plot(denver.data2.1, "obs", ID="d_43", xlab="", ylab="BC (log ug/m3)",
     main="Temporal trend d_43")
plot(denver.data2.1, "res", ID="d_43", xlab="", ylab="BC (log ug/m3)")
plot(denver.data2.1, "acf", ID="d_43")
plot(denver.data2.1, "pacf", ID="d_43")

par(mar=c(3.3,3.3,1.5,1), mgp=c(2,1,0))
layout(matrix(c(1,1,2,2,3,4), 3, 2, byrow=TRUE))
plot(denver.data2.1, "obs", ID="d_53", xlab="", ylab="BC (log ug/m3)",
     main="Temporal trend d_53")
plot(denver.data2.1, "res", ID="d_53", xlab="", ylab="BC (log ug/m3)")
plot(denver.data2.1, "acf", ID="d_53")
plot(denver.data2.1, "pacf", ID="d_53")

par(mar=c(3.3,3.3,1.5,1), mgp=c(2,1,0))
layout(matrix(c(1,1,2,2,3,4), 3, 2, byrow=TRUE))
plot(denver.data2.1, "obs", ID="central", xlab="", ylab="BC (log ug/m3)",
     main="Temporal trend central")
plot(denver.data2.1, "res", ID="central", xlab="", ylab="BC (log ug/m3)")
plot(denver.data2.1, "acf", ID="central")
plot(denver.data2.1, "pacf", ID="central")

2.5.2 1 basis function, 8 df per year

denver.data2.1a <- denver.data
denver.data2.1a$trend <- pol_STdata1a$trend
print(denver.data2.1a)
## STdata-object with:
##  No. locations: 69 (observed: 64)
##  No. time points: 622 (observed: 231)
##  No. obs: 1021
## 
## Trend with 1 basis function(s):
## [1] "V1"
## with dates:
##  2008-12-29 to 2020-12-07
## 
## 19 covariate(s):
##  [1] "ID"              "lon"             "lat"             "tree_cover_100" 
##  [5] "impervious_2500" "low_int_100"     "med_int_50"      "med_int_2500"   
##  [9] "high_int_100"    "pop_den_50"      "dist_m_cafos"    "dist_m_og_wells"
## [13] "dist_m_wwtp"     "aadt_100"        "aadt_250"        "aadt_1000"      
## [17] "x"               "y"               "type"           
## 
## 2 spatio-temporal covariate(s):
## [1] "bc_st_no2" "bc_st_smk"
## 
## All sites:
## central    dist 
##       1      68 
## Observed:
## central    dist 
##       1      63 
## 
## For central:
##   Number of obs: 231
##   Dates: 2016-02-08 to 2020-11-30
## For dist:
##   Number of obs: 790
##   Dates: 2018-05-07 to 2020-04-06
par(mar=c(3.3,3.3,1.5,1), mgp=c(2,1,0))
layout(matrix(c(1,1,2,2,3,4), 3, 2, byrow=TRUE))
plot(denver.data2.1a, "obs", ID="d_1", xlab="", ylab="BC (log ug/m3)",
     main="Temporal trend d_1")
plot(denver.data2.1a, "res", ID="d_1", xlab="", ylab="BC (log ug/m3)")
plot(denver.data2.1a, "acf", ID="d_1")
plot(denver.data2.1a, "pacf", ID="d_1")

par(mar=c(3.3,3.3,1.5,1), mgp=c(2,1,0))
layout(matrix(c(1,1,2,2,3,4), 3, 2, byrow=TRUE))
plot(denver.data2.1a, "obs", ID="d_10", xlab="", ylab="BC (log ug/m3)",
     main="Temporal trend d_43")
plot(denver.data2.1a, "res", ID="d_10", xlab="", ylab="BC (log ug/m3)")
plot(denver.data2.1a, "acf", ID="d_10")
plot(denver.data2.1a, "pacf", ID="d_10")

par(mar=c(3.3,3.3,1.5,1), mgp=c(2,1,0))
layout(matrix(c(1,1,2,2,3,4), 3, 2, byrow=TRUE))
plot(denver.data2.1a, "obs", ID="d_43", xlab="", ylab="BC (log ug/m3)",
     main="Temporal trend d_43")
plot(denver.data2.1a, "res", ID="d_43", xlab="", ylab="BC (log ug/m3)")
plot(denver.data2.1a, "acf", ID="d_43")
plot(denver.data2.1a, "pacf", ID="d_43")

par(mar=c(3.3,3.3,1.5,1), mgp=c(2,1,0))
layout(matrix(c(1,1,2,2,3,4), 3, 2, byrow=TRUE))
plot(denver.data2.1a, "obs", ID="d_53", xlab="", ylab="BC (log ug/m3)",
     main="Temporal trend d_53")
plot(denver.data2.1a, "res", ID="d_53", xlab="", ylab="BC (log ug/m3)")
plot(denver.data2.1a, "acf", ID="d_53")
plot(denver.data2.1a, "pacf", ID="d_53")

par(mar=c(3.3,3.3,1.5,1), mgp=c(2,1,0))
layout(matrix(c(1,1,2,2,3,4), 3, 2, byrow=TRUE))
plot(denver.data2.1a, "obs", ID="central", xlab="", ylab="BC (log ug/m3)",
     main="Temporal trend central")
plot(denver.data2.1a, "res", ID="central", xlab="", ylab="BC (log ug/m3)")
plot(denver.data2.1a, "acf", ID="central")
plot(denver.data2.1a, "pacf", ID="central")

2.5.3 2 basis functions, 4 df per year

denver.data2.2 <- denver.data
denver.data2.2$trend <- pol_STdata2$trend
print(denver.data2.2)
## STdata-object with:
##  No. locations: 69 (observed: 64)
##  No. time points: 622 (observed: 231)
##  No. obs: 1021
## 
## Trend with 2 basis function(s):
## [1] "V1" "V2"
## with dates:
##  2008-12-29 to 2020-12-07
## 
## 19 covariate(s):
##  [1] "ID"              "lon"             "lat"             "tree_cover_100" 
##  [5] "impervious_2500" "low_int_100"     "med_int_50"      "med_int_2500"   
##  [9] "high_int_100"    "pop_den_50"      "dist_m_cafos"    "dist_m_og_wells"
## [13] "dist_m_wwtp"     "aadt_100"        "aadt_250"        "aadt_1000"      
## [17] "x"               "y"               "type"           
## 
## 2 spatio-temporal covariate(s):
## [1] "bc_st_no2" "bc_st_smk"
## 
## All sites:
## central    dist 
##       1      68 
## Observed:
## central    dist 
##       1      63 
## 
## For central:
##   Number of obs: 231
##   Dates: 2016-02-08 to 2020-11-30
## For dist:
##   Number of obs: 790
##   Dates: 2018-05-07 to 2020-04-06
par(mar=c(3.3,3.3,1.5,1), mgp=c(2,1,0))
layout(matrix(c(1,1,2,2,3,4), 3, 2, byrow=TRUE))
plot(denver.data2.2, "obs", ID="d_1", xlab="", ylab="BC (log ug/m3)",
     main="Temporal trend d_1")
plot(denver.data2.2, "res", ID="d_1", xlab="", ylab="BC (log ug/m3)")
plot(denver.data2.2, "acf", ID="d_1")
plot(denver.data2.2, "pacf", ID="d_1")

par(mar=c(3.3,3.3,1.5,1), mgp=c(2,1,0))
layout(matrix(c(1,1,2,2,3,4), 3, 2, byrow=TRUE))
plot(denver.data2.2, "obs", ID="d_10", xlab="", ylab="BC (log ug/m3)",
     main="Temporal trend d_10")
plot(denver.data2.2, "res", ID="d_10", xlab="", ylab="BC (log ug/m3)")
plot(denver.data2.2, "acf", ID="d_10")
plot(denver.data2.2, "pacf", ID="d_10")

par(mar=c(3.3,3.3,1.5,1), mgp=c(2,1,0))
layout(matrix(c(1,1,2,2,3,4), 3, 2, byrow=TRUE))
plot(denver.data2.2, "obs", ID="d_43", xlab="", ylab="BC (log ug/m3)",
     main="Temporal trend d_43")
plot(denver.data2.2, "res", ID="d_43", xlab="", ylab="BC (log ug/m3)")
plot(denver.data2.2, "acf", ID="d_43")
plot(denver.data2.2, "pacf", ID="d_43")

par(mar=c(3.3,3.3,1.5,1), mgp=c(2,1,0))
layout(matrix(c(1,1,2,2,3,4), 3, 2, byrow=TRUE))
plot(denver.data2.2, "obs", ID="d_53", xlab="", ylab="BC (log ug/m3)",
     main="Temporal trend d_53")
plot(denver.data2.2, "res", ID="d_53", xlab="", ylab="BC (log ug/m3)")
plot(denver.data2.2, "acf", ID="d_53")
plot(denver.data2.2, "pacf", ID="d_53")

par(mar=c(3.3,3.3,1.5,1), mgp=c(2,1,0))
layout(matrix(c(1,1,2,2,3,4), 3, 2, byrow=TRUE))
plot(denver.data2.2, "obs", ID="central", xlab="", ylab="BC (log ug/m3)",
     main="Temporal trend central")
plot(denver.data2.2, "res", ID="central", xlab="", ylab="BC (log ug/m3)")
plot(denver.data2.2, "acf", ID="central")
plot(denver.data2.2, "pacf", ID="central")

2.5.4 2 basis functions, 8 df per year

denver.data2.2a <- denver.data
denver.data2.2a$trend <- pol_STdata2a$trend
print(denver.data2.2a)
## STdata-object with:
##  No. locations: 69 (observed: 64)
##  No. time points: 622 (observed: 231)
##  No. obs: 1021
## 
## Trend with 2 basis function(s):
## [1] "V1" "V2"
## with dates:
##  2008-12-29 to 2020-12-07
## 
## 19 covariate(s):
##  [1] "ID"              "lon"             "lat"             "tree_cover_100" 
##  [5] "impervious_2500" "low_int_100"     "med_int_50"      "med_int_2500"   
##  [9] "high_int_100"    "pop_den_50"      "dist_m_cafos"    "dist_m_og_wells"
## [13] "dist_m_wwtp"     "aadt_100"        "aadt_250"        "aadt_1000"      
## [17] "x"               "y"               "type"           
## 
## 2 spatio-temporal covariate(s):
## [1] "bc_st_no2" "bc_st_smk"
## 
## All sites:
## central    dist 
##       1      68 
## Observed:
## central    dist 
##       1      63 
## 
## For central:
##   Number of obs: 231
##   Dates: 2016-02-08 to 2020-11-30
## For dist:
##   Number of obs: 790
##   Dates: 2018-05-07 to 2020-04-06
par(mar=c(3.3,3.3,1.5,1), mgp=c(2,1,0))
layout(matrix(c(1,1,2,2,3,4), 3, 2, byrow=TRUE))
plot(denver.data2.2a, "obs", ID="d_1", xlab="", ylab="BC (log ug/m3)",
     main="Temporal trend d_1")
plot(denver.data2.2a, "res", ID="d_1", xlab="", ylab="BC (log ug/m3)")
plot(denver.data2.2a, "acf", ID="d_1")
plot(denver.data2.2a, "pacf", ID="d_1")

par(mar=c(3.3,3.3,1.5,1), mgp=c(2,1,0))
layout(matrix(c(1,1,2,2,3,4), 3, 2, byrow=TRUE))
plot(denver.data2.2a, "obs", ID="d_10", xlab="", ylab="BC (log ug/m3)",
     main="Temporal trend d_10")
plot(denver.data2.2a, "res", ID="d_10", xlab="", ylab="BC (log ug/m3)")
plot(denver.data2.2a, "acf", ID="d_10")
plot(denver.data2.2a, "pacf", ID="d_10")

par(mar=c(3.3,3.3,1.5,1), mgp=c(2,1,0))
layout(matrix(c(1,1,2,2,3,4), 3, 2, byrow=TRUE))
plot(denver.data2.2a, "obs", ID="d_43", xlab="", ylab="BC (log ug/m3)",
     main="Temporal trend d_43")
plot(denver.data2.2a, "res", ID="d_43", xlab="", ylab="BC (log ug/m3)")
plot(denver.data2.2a, "acf", ID="d_43")
plot(denver.data2.2a, "pacf", ID="d_43")

par(mar=c(3.3,3.3,1.5,1), mgp=c(2,1,0))
layout(matrix(c(1,1,2,2,3,4), 3, 2, byrow=TRUE))
plot(denver.data2.2a, "obs", ID="d_53", xlab="", ylab="BC (log ug/m3)",
     main="Temporal trend d_53")
plot(denver.data2.2a, "res", ID="d_53", xlab="", ylab="BC (log ug/m3)")
plot(denver.data2.2a, "acf", ID="d_53")
plot(denver.data2.2a, "pacf", ID="d_53")

par(mar=c(3.3,3.3,1.5,1), mgp=c(2,1,0))
layout(matrix(c(1,1,2,2,3,4), 3, 2, byrow=TRUE))
plot(denver.data2.2a, "obs", ID="central", xlab="", ylab="BC (log ug/m3)",
     main="Temporal trend central")
plot(denver.data2.2a, "res", ID="central", xlab="", ylab="BC (log ug/m3)")
plot(denver.data2.2a, "acf", ID="central")
plot(denver.data2.2a, "pacf", ID="central")

3 Testing possible models

Building on the results of the previous modeling efforts, I’m going to update “Model 4.2”, which used an ’exp covariance structure for the error term and iid for the beta fields.

3.1 Model A: 1 temporal trend (basis) function (4 df per year)

3.1.1 Create the model object

For this version of the model, use iid for cov.beta (beta, beta1) and exp for cov.nu (error). Here we can specify different LUR formluae. The length of the LUR list should be number of basis functions + 1.

denver.data.A <- denver.data2.1
names(denver.data.A$covars)
##  [1] "ID"              "lon"             "lat"             "tree_cover_100" 
##  [5] "impervious_2500" "low_int_100"     "med_int_50"      "med_int_2500"   
##  [9] "high_int_100"    "pop_den_50"      "dist_m_cafos"    "dist_m_og_wells"
## [13] "dist_m_wwtp"     "aadt_100"        "aadt_250"        "aadt_1000"      
## [17] "x"               "y"               "type"
LUR.A <- list(covar_fun, covar_fun)

cov.beta.A <-  list(covf = c("iid", "iid"), nugget = T)
cov.nu.A <- list(covf = "exp", nugget = T, random.effect = FALSE)
locations.A <- list(coords = c("x","y"), long.lat = c("lon","lat"), others= "type")

denver.model.A <- createSTmodel(denver.data.A, LUR = LUR.A,
                                ST = c("bc_st_no2", "bc_st_smk"),
                                cov.beta = cov.beta.A, cov.nu = cov.nu.A,
                                locations = locations.A)
## No trend $trend.fnc object detected, STdata probably from old version of the package.
## $trend.fnc has been added based on spline fit to elements in STmodel$trend.
denver.model.A
## STmodel-object with:
##  No. locations: 69 (observed: 64)
##  No. time points: 622 (observed: 231)
##  No. obs: 1021
## 
## Trend with 1 basis function(s):
## [1] "V1"
## with dates:
##  2008-12-29 to 2020-12-07
## 
## Models for the beta-fields are:
## $const
## ~~tree_cover_100 + impervious_2500 + low_int_100 + med_int_50 + 
##     med_int_2500 + high_int_100 + pop_den_50 + dist_m_cafos + 
##     dist_m_og_wells + dist_m_wwtp + aadt_100 + aadt_250 + aadt_1000
## 
## $V1
## ~~tree_cover_100 + impervious_2500 + low_int_100 + med_int_50 + 
##     med_int_2500 + high_int_100 + pop_den_50 + dist_m_cafos + 
##     dist_m_og_wells + dist_m_wwtp + aadt_100 + aadt_250 + aadt_1000
## 
## 2 spatio-temporal covariate(s):
## [1] "bc_st_no2" "bc_st_smk"
## 
## Covariance model for the beta-field(s):
##  Covariance type(s): iid, iid 
##  Nugget: Yes, Yes 
## Covariance model for the nu-field(s):
##  Covariance type: exp 
##  Nugget: ~1 
##  Random effect: No 
## All sites:
## central    dist 
##       1      68 
## Observed:
## central    dist 
##       1      63 
## 
## For central:
##   Number of obs: 231
##   Dates: 2016-02-08 to 2020-11-30
## For dist:
##   Number of obs: 790
##   Dates: 2018-05-07 to 2020-04-06

3.1.2 Estimate model parameters

  • nugget values can range from -1 to -6
  • range values can range from 0 to 4
  • try starting values where the sill and nugget are the same (e.g., both 0) and where they are very different (e.g., 0 and -5)
  • make sure the starting values are pretty different
set.seed(1000)

names <- loglikeSTnames(denver.model.A, all=FALSE)
names
## [1] "log.nugget.const.iid"          "log.nugget.V1.iid"            
## [3] "nu.log.range.exp"              "nu.log.sill.exp"              
## [5] "nu.log.nugget.(Intercept).exp"
# x.init.A <- cbind(c(0, 0, 0, 0, 0),
#                   c(-1, -1, 0, -1, -1),
#                   c(-1, -1, 0, -5, -1),
#                   c(-5, -5, 0, -1, -5),
#                   c(-5, -5, 0, -5, -5),
#                   c(-1, -1, 2, -1, -1),
#                   c(-1, -1, 2, -5, -1),
#                   c(-5, -5, 2, -1, -5),
#                   c(-5, -5, 2, -5, -5),
#                   c(-1, -1, 4, -1, -1),
#                   c(-1, -1, 4, -5, -1),
#                   c(-5, -5, 4, -1, -5),
#                   c(-5, -5, 4, -5, -5))

x.init.A <- cbind(c(0, 0, 0, 0, 0),
                  c(-1, -1, 4, -5, -1),
                  c(-5, -5, 4, -1, -5),
                  c(-1, -1, 8, -5, -1),
                  c(-5, -5, 8, -1, -5),
                  c(-1, -1, 12, -5, -1),
                  c(-5, -5, 12, -1, -5),
                  c(-7, -7, 12, -3, -3))

rownames(x.init.A) <- loglikeSTnames(denver.model.A, all=FALSE)
x.init.A
##                               [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
## log.nugget.const.iid             0   -1   -5   -1   -5   -1   -5   -7
## log.nugget.V1.iid                0   -1   -5   -1   -5   -1   -5   -7
## nu.log.range.exp                 0    4    4    8    8   12   12   12
## nu.log.sill.exp                  0   -5   -1   -5   -1   -5   -1   -3
## nu.log.nugget.(Intercept).exp    0   -1   -5   -1   -5   -1   -5   -3
#' Difference from tutorial: use Josh Keller's version of the function
source(here::here("Code", "functions_model.R"))
est.denver.model.A <- estimate.STmodel(denver.model.A, x.init.A)
## Optimisation using starting value 1/8
## N = 5, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate     0  f=       459.08  |proj g|=           15
## At iterate    10  f =      -1309.8  |proj g|=        1.2358
## At iterate    20  f =        -1310  |proj g|=     0.0064564
## 
## iterations 21
## function evaluations 35
## segments explored during Cauchy searches 25
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 2
## norm of the final projected gradient 0.00646143
## final function value -1310.01
## 
## F = -1310.01
## final  value -1310.007198 
## converged
## Optimisation using starting value 2/8
## N = 5, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate     0  f=      -331.87  |proj g|=           14
## At iterate    10  f =      -1306.2  |proj g|=         1.925
## ys=-2.142e-01  -gs= 4.291e-02, BFGS update SKIPPED
## At iterate    20  f =      -1307.4  |proj g|=        11.314
## At iterate    30  f =        -1310  |proj g|=        0.5797
## At iterate    40  f =      -1402.8  |proj g|=        20.928
## At iterate    50  f =      -1572.4  |proj g|=      0.014546
## 
## iterations 51
## function evaluations 93
## segments explored during Cauchy searches 54
## BFGS updates skipped 1
## active bounds at final generalized Cauchy point 0
## norm of the final projected gradient 0.014546
## final function value -1572.4
## 
## F = -1572.4
## final  value -1572.396317 
## converged
## Optimisation using starting value 3/8
## N = 5, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate     0  f=      -458.66  |proj g|=           14
## At iterate    10  f =      -1304.7  |proj g|=        1.9297
## At iterate    20  f =      -1308.9  |proj g|=         10.68
## ys=-2.936e+00  -gs= 9.218e-01, BFGS update SKIPPED
## At iterate    30  f =        -1540  |proj g|=        20.292
## 
## iterations 39
## function evaluations 53
## segments explored during Cauchy searches 40
## BFGS updates skipped 1
## active bounds at final generalized Cauchy point 0
## norm of the final projected gradient 0.00692961
## final function value -1572.4
## 
## F = -1572.4
## final  value -1572.397522 
## converged
## Optimisation using starting value 4/8
## N = 5, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate     0  f=      -332.17  |proj g|=           14
## At iterate    10  f =      -1571.2  |proj g|=        3.1726
## At iterate    20  f =      -1576.3  |proj g|=       0.32053
## Bad direction in the line search;
##    refresh the lbfgs memory and restart the iteration.
## 
## iterations 29
## function evaluations 79
## segments explored during Cauchy searches 33
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 0
## norm of the final projected gradient 0.00351314
## final function value -1576.3
## 
## F = -1576.3
## Warning:  more than 10 function and gradient evaluations
##    in the last line search
## final  value -1576.303657 
## converged
## Optimisation using starting value 5/8
## N = 5, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate     0  f=      -546.73  |proj g|=           14
## At iterate    10  f =        -1576  |proj g|=        1.1049
## At iterate    20  f =      -1576.3  |proj g|=       0.08239
## 
## iterations 24
## function evaluations 46
## segments explored during Cauchy searches 26
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 0
## norm of the final projected gradient 0.0133713
## final function value -1576.3
## 
## F = -1576.3
## Warning:  more than 10 function and gradient evaluations
##    in the last line search
## final  value -1576.303693 
## converged
## Optimisation using starting value 6/8
## N = 5, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate     0  f=      -335.65  |proj g|=           14
## At iterate    10  f =      -1576.1  |proj g|=        6.8273
## 
## iterations 19
## function evaluations 37
## segments explored during Cauchy searches 22
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 0
## norm of the final projected gradient 0.076101
## final function value -1576.3
## 
## F = -1576.3
## Warning:  more than 10 function and gradient evaluations
##    in the last line search
## final  value -1576.303485 
## converged
## Optimisation using starting value 7/8
## N = 5, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate     0  f=      -1296.8  |proj g|=           14
## At iterate    10  f =      -1576.3  |proj g|=       0.36275
## At iterate    20  f =      -1576.3  |proj g|=       0.06807
## 
## iterations 26
## function evaluations 35
## segments explored during Cauchy searches 30
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 0
## norm of the final projected gradient 0.0533678
## final function value -1576.3
## 
## F = -1576.3
## final  value -1576.303711 
## converged
## Optimisation using starting value 8/8
## N = 5, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate     0  f=      -1228.8  |proj g|=           12
## At iterate    10  f =      -1576.3  |proj g|=        1.8747
## At iterate    20  f =      -1576.3  |proj g|=      0.082774
## 
## iterations 22
## function evaluations 28
## segments explored during Cauchy searches 25
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 0
## norm of the final projected gradient 0.0441012
## final function value -1576.3
## 
## F = -1576.3
## final  value -1576.303691 
## converged
print(est.denver.model.A)
## Optimisation for STmodel with 8 starting points.
##   Results: 5 converged, 3 not converged, 0 failed.
##   Best result for starting point 7, optimisation has converged
## 
## No fixed parameters.
## 
## Estimated parameters for all starting point(s):
##                                           [,1]          [,2]          [,3]
## gamma.bc_st_no2                 0.018662027683   0.011947334   0.011948431
## gamma.bc_st_smk                -0.045770923734  -0.024770380  -0.024768403
## alpha.const.(Intercept)        -0.144213749735   0.043061298   0.043040010
## alpha.const.tree_cover_100     -0.001271761777   0.003719008   0.003718359
## alpha.const.impervious_2500     0.052029988535   0.044170288   0.044168348
## alpha.const.low_int_100        -0.001653106787   0.009016166   0.009015446
## alpha.const.med_int_50         -0.002288302995   0.002734351   0.002734279
## alpha.const.med_int_2500       -0.010397390358  -0.017786781  -0.017786385
## alpha.const.high_int_100        0.010312389425   0.016163245   0.016162803
## alpha.const.pop_den_50         -0.015664209451  -0.013249908  -0.013248458
## alpha.const.dist_m_cafos       -0.045208992335  -0.037819717  -0.037820572
## alpha.const.dist_m_og_wells     0.010367781069   0.002794518   0.002795706
## alpha.const.dist_m_wwtp         0.005191651255   0.006739532   0.006738714
## alpha.const.aadt_100            0.034921854611   0.034342561   0.034342776
## alpha.const.aadt_250            0.008441674556   0.006778407   0.006777233
## alpha.const.aadt_1000          -0.001209607701   0.001401806   0.001401659
## alpha.V1.(Intercept)            0.186004760519   0.222990414   0.222990813
## alpha.V1.tree_cover_100        -0.020164715529  -0.016421057  -0.016422865
## alpha.V1.impervious_2500       -0.020298784866  -0.038355053  -0.038358586
## alpha.V1.low_int_100           -0.016748233869  -0.008141071  -0.008142264
## alpha.V1.med_int_50            -0.022951471286  -0.017767597  -0.017767805
## alpha.V1.med_int_2500           0.026917676047   0.025099230   0.025101869
## alpha.V1.high_int_100           0.004004735689   0.009705730   0.009702999
## alpha.V1.pop_den_50            -0.000008909438   0.002641688   0.002640935
## alpha.V1.dist_m_cafos           0.011071875338   0.007011703   0.007010877
## alpha.V1.dist_m_og_wells       -0.004784189400   0.001697012   0.001696335
## alpha.V1.dist_m_wwtp           -0.009890449702  -0.011769544  -0.011770430
## alpha.V1.aadt_100              -0.046053857986  -0.041872920  -0.041868833
## alpha.V1.aadt_250               0.009842151159   0.007113844   0.007113932
## alpha.V1.aadt_1000             -0.007732179733  -0.005245424  -0.005245087
## log.nugget.const.iid          -15.000000000000 -14.999371995 -14.998712098
## log.nugget.V1.iid             -15.000000000000 -14.999810884 -14.463366983
## nu.log.range.exp                0.000000000000  12.464310836  12.463721125
## nu.log.sill.exp                -4.390684341790  -3.022719588  -3.022834153
## nu.log.nugget.(Intercept).exp  -4.135497997773  -4.771806577  -4.771972935
##                                       [,4]         [,5]         [,6]
## gamma.bc_st_no2                0.012195145  0.012191546  0.012201358
## gamma.bc_st_smk               -0.023368096 -0.023376519 -0.023355592
## alpha.const.(Intercept)        0.034510427  0.034586281  0.034389759
## alpha.const.tree_cover_100     0.002813683  0.002817342  0.002809626
## alpha.const.impervious_2500    0.044148457  0.044155063  0.044138336
## alpha.const.low_int_100        0.008272737  0.008276440  0.008268586
## alpha.const.med_int_50         0.001462357  0.001462844  0.001465164
## alpha.const.med_int_2500      -0.019261212 -0.019261589 -0.019257039
## alpha.const.high_int_100       0.014684016  0.014687448  0.014683691
## alpha.const.pop_den_50        -0.011372682 -0.011379302 -0.011366384
## alpha.const.dist_m_cafos      -0.038056720 -0.038053453 -0.038064354
## alpha.const.dist_m_og_wells    0.002481102  0.002476676  0.002492421
## alpha.const.dist_m_wwtp        0.005861502  0.005865456  0.005856464
## alpha.const.aadt_100           0.033664281  0.033663797  0.033666382
## alpha.const.aadt_250           0.007104613  0.007108579  0.007095515
## alpha.const.aadt_1000          0.002106792  0.002106273  0.002102821
## alpha.V1.(Intercept)           0.226073249  0.226100754  0.226019975
## alpha.V1.tree_cover_100       -0.016992745 -0.016986819 -0.017002648
## alpha.V1.impervious_2500      -0.036778794 -0.036765912 -0.036806336
## alpha.V1.low_int_100          -0.008874293 -0.008869243 -0.008882278
## alpha.V1.med_int_50           -0.016196431 -0.016196237 -0.016203423
## alpha.V1.med_int_2500          0.024306189  0.024301028  0.024314867
## alpha.V1.high_int_100          0.009012393  0.009021579  0.008997412
## alpha.V1.pop_den_50            0.002208200  0.002208468  0.002210902
## alpha.V1.dist_m_cafos          0.006169577  0.006170451  0.006173326
## alpha.V1.dist_m_og_wells       0.001844731  0.001849461  0.001832932
## alpha.V1.dist_m_wwtp          -0.010408832 -0.010408312 -0.010415084
## alpha.V1.aadt_100             -0.040115551 -0.040129281 -0.040090907
## alpha.V1.aadt_250              0.007671507  0.007668694  0.007672913
## alpha.V1.aadt_1000            -0.005173673 -0.005174620 -0.005171007
## log.nugget.const.iid          -7.797159485 -7.798415070 -7.800270351
## log.nugget.V1.iid             -7.521046376 -7.525829277 -7.515414412
## nu.log.range.exp              12.444970355 12.446646679 12.441842084
## nu.log.sill.exp               -3.008558413 -3.008570123 -3.008497539
## nu.log.nugget.(Intercept).exp -4.861112033 -4.860671802 -4.861403855
##                                       [,7]         [,8]
## gamma.bc_st_no2                0.012187315  0.012191285
## gamma.bc_st_smk               -0.023383491 -0.023379810
## alpha.const.(Intercept)        0.034674416  0.034601251
## alpha.const.tree_cover_100     0.002819106  0.002814760
## alpha.const.impervious_2500    0.044164671  0.044151847
## alpha.const.low_int_100        0.008279085  0.008274175
## alpha.const.med_int_50         0.001460604  0.001462762
## alpha.const.med_int_2500      -0.019266096 -0.019261523
## alpha.const.high_int_100       0.014687926  0.014685727
## alpha.const.pop_den_50        -0.011382275 -0.011375796
## alpha.const.dist_m_cafos      -0.038050220 -0.038055208
## alpha.const.dist_m_og_wells    0.002471527  0.002480366
## alpha.const.dist_m_wwtp        0.005867949  0.005862805
## alpha.const.aadt_100           0.033662812  0.033665304
## alpha.const.aadt_250           0.007114408  0.007105814
## alpha.const.aadt_1000          0.002106229  0.002104928
## alpha.V1.(Intercept)           0.226138646  0.226089325
## alpha.V1.tree_cover_100       -0.016979945 -0.016990544
## alpha.V1.impervious_2500      -0.036741674 -0.036773466
## alpha.V1.low_int_100          -0.008864162 -0.008872411
## alpha.V1.med_int_50           -0.016194695 -0.016197912
## alpha.V1.med_int_2500          0.024284973  0.024299875
## alpha.V1.high_int_100          0.009035275  0.009017501
## alpha.V1.pop_den_50            0.002211442  0.002210976
## alpha.V1.dist_m_cafos          0.006174037  0.006172359
## alpha.V1.dist_m_og_wells       0.001851383  0.001843365
## alpha.V1.dist_m_wwtp          -0.010402080 -0.010408105
## alpha.V1.aadt_100             -0.040142669 -0.040117310
## alpha.V1.aadt_250              0.007667448  0.007669650
## alpha.V1.aadt_1000            -0.005177372 -0.005173688
## log.nugget.const.iid          -7.797036706 -7.798340467
## log.nugget.V1.iid             -7.526945504 -7.521553058
## nu.log.range.exp              12.449157106 12.446308860
## nu.log.sill.exp               -3.008551080 -3.007960669
## nu.log.nugget.(Intercept).exp -4.860112385 -4.861008888
## 
## Function value(s):
## [1] 1310.007 1572.396 1572.398 1576.304 1576.304 1576.303 1576.304 1576.304

3.1.3 Cross-validation

Define the CV groups

set.seed(1000)
site_idsA <- colnames(denver.model.A$D.beta)[which(str_detect(colnames(denver.model.A$D.beta), "d_") == T)]

unique(colnames(bc_obs2))
##  [1] "central" "d_1"     "d_10"    "d_11"    "d_12"    "d_13"    "d_14"   
##  [8] "d_15"    "d_16"    "d_17"    "d_18"    "d_19"    "d_2"     "d_20"   
## [15] "d_21"    "d_22"    "d_23"    "d_24"    "d_25"    "d_26"    "d_27"   
## [22] "d_28"    "d_29"    "d_3"     "d_30"    "d_31"    "d_32"    "d_33"   
## [29] "d_34"    "d_35"    "d_36"    "d_37"    "d_38"    "d_39"    "d_4"    
## [36] "d_40"    "d_41"    "d_42"    "d_43"    "d_44"    "d_45"    "d_46"   
## [43] "d_47"    "d_48"    "d_49"    "d_5"     "d_50"    "d_51"    "d_52"   
## [50] "d_53"    "d_54"    "d_55"    "d_56"    "d_57"    "d_58"    "d_59"   
## [57] "d_6"     "d_60"    "d_61"    "d_62"    "d_63"    "d_64"    "d_65"   
## [64] "d_66"    "d_67"    "d_68"    "d_7"     "d_8"     "d_9"
Ind.cv.A <- createCV(denver.model.A, groups = 10, min.dist = .1,
                     subset = site_idsA)

ID.cv.A <- sapply(split(denver.model.A$obs$ID, Ind.cv.A), unique)
print(sapply(ID.cv.A, length))
##  0  1  2  3  4  5  6  7  8  9 10 
##  1  7  7  7  6  6  6  6  6  6  6
table(Ind.cv.A)
## Ind.cv.A
##   0   1   2   3   4   5   6   7   8   9  10 
## 231  75  80  97  73  64  93  80  81  66  81
I.col.A <- apply(sapply(ID.cv.A, function(x) denver.model.A$locations$ID%in% x), 1,
                       function(x) if(sum(x)==1) which(x) else 0)
names(I.col.A) <- denver.model.A$locations$ID
print(I.col.A)
## central     d_1    d_10    d_11    d_12    d_13    d_14    d_15    d_16    d_17 
##       1       6       2       3       4       4       4       9       4      10 
##    d_18    d_19     d_2    d_20    d_21    d_22    d_23    d_25    d_26    d_28 
##       8       5      11       3       2      11       4       9      11       9 
##    d_29     d_3    d_31    d_32    d_33    d_34    d_35    d_36    d_37    d_38 
##       2       4       6       8       5       5      10       8       3       8 
##    d_39     d_4    d_40    d_41    d_42    d_43    d_44    d_45    d_46    d_47 
##       6       7       7       3       7       8      11       3       9       5 
##    d_48    d_49     d_5    d_50    d_51    d_52    d_53    d_54    d_55    d_56 
##       9       4       6       9      10       3       7       5       7       6 
##    d_57    d_58     d_6    d_60    d_61    d_62    d_63    d_64    d_65    d_66 
##       2       2      11       5       2      11      10       3       8       6 
##    d_67    d_68     d_7     d_8    d_24    d_27    d_30    d_59     d_9 
##       7      10       2      10       0       0       0       0       0
par(mfrow=c(1,1))
plot(denver.model.A$locations$long,
     denver.model.A$locations$lat,
     pch=23+floor(I.col.A/max(I.col.A)+.5), bg=I.col.A,
     xlab="Longitude", ylab="Latitude")

map("county", "colorado", col="#FFFF0055",fill=TRUE, add=TRUE)
## [[1]]
## NULL

ID starting conditions using previous model without CV:

x.init.A.cv <- coef(est.denver.model.A, pars="cov")[,c("par","init")]
x.init.A.cv

Run the model with cross validation

est.denver.A.cv <- estimateCV(denver.model.A, x.init.A.cv, Ind.cv.A)
## 
## ***************************
## Estimation of CV-set 1/10
## Optimisation using starting value 1/2
## N = 5, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate     0  f=      -1430.9  |proj g|=       11.381
## At iterate    10  f =      -1431.7  |proj g|=       0.27352
## 
## iterations 17
## function evaluations 30
## segments explored during Cauchy searches 17
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 0
## norm of the final projected gradient 0.0379531
## final function value -1431.71
## 
## F = -1431.71
## final  value -1431.711037 
## converged
## 
## Optimisation using starting value 2/2
## N = 5, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate     0  f=      -1163.8  |proj g|=           14
## At iterate    10  f =      -1431.5  |proj g|=       0.61862
## At iterate    20  f =      -1431.7  |proj g|=       0.13907
## Bad direction in the line search;
##    refresh the lbfgs memory and restart the iteration.
## 
## iterations 25
## function evaluations 61
## segments explored during Cauchy searches 30
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 0
## norm of the final projected gradient 0.0694925
## final function value -1431.71
## 
## F = -1431.71
## Warning:  more than 10 function and gradient evaluations
##    in the last line search
## final  value -1431.711043 
## converged
## 
## 
## ***************************
## Estimation of CV-set 2/10
## Optimisation using starting value 1/2
## N = 5, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate     0  f=      -1435.9  |proj g|=       1.0766
## At iterate    10  f =      -1435.9  |proj g|=      0.020506
## 
## iterations 10
## function evaluations 28
## segments explored during Cauchy searches 10
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 0
## norm of the final projected gradient 0.0205057
## final function value -1435.92
## 
## F = -1435.92
## Warning:  more than 10 function and gradient evaluations
##    in the last line search
## final  value -1435.924803 
## converged
## 
## Optimisation using starting value 2/2
## N = 5, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate     0  f=      -1169.1  |proj g|=           14
## At iterate    10  f =      -1435.9  |proj g|=       0.22268
## 
## iterations 19
## function evaluations 34
## segments explored during Cauchy searches 23
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 0
## norm of the final projected gradient 0.00810512
## final function value -1435.92
## 
## F = -1435.92
## final  value -1435.924810 
## converged
## 
## 
## ***************************
## Estimation of CV-set 3/10
## Optimisation using starting value 1/2
## N = 5, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate     0  f=      -1419.9  |proj g|=        10.14
## At iterate    10  f =      -1420.6  |proj g|=       0.23026
## 
## iterations 17
## function evaluations 31
## segments explored during Cauchy searches 19
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 0
## norm of the final projected gradient 0.0225161
## final function value -1420.64
## 
## F = -1420.64
## final  value -1420.637129 
## converged
## 
## Optimisation using starting value 2/2
## N = 5, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate     0  f=      -1149.9  |proj g|=           14
## At iterate    10  f =      -1420.6  |proj g|=       0.40793
## At iterate    20  f =      -1420.6  |proj g|=       0.03382
## 
## iterations 20
## function evaluations 25
## segments explored during Cauchy searches 24
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 0
## norm of the final projected gradient 0.0338201
## final function value -1420.64
## 
## F = -1420.64
## final  value -1420.637053 
## converged
## 
## 
## ***************************
## Estimation of CV-set 4/10
## Optimisation using starting value 1/2
## N = 5, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate     0  f=      -1425.1  |proj g|=       17.198
## At iterate    10  f =      -1425.9  |proj g|=       0.20796
## 
## iterations 17
## function evaluations 30
## segments explored during Cauchy searches 17
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 0
## norm of the final projected gradient 0.00811343
## final function value -1425.92
## 
## F = -1425.92
## final  value -1425.924348 
## converged
## 
## Optimisation using starting value 2/2
## N = 5, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate     0  f=      -1161.3  |proj g|=           14
## At iterate    10  f =      -1425.7  |proj g|=       0.64768
## At iterate    20  f =      -1425.9  |proj g|=     0.0094272
## Bad direction in the line search;
##    refresh the lbfgs memory and restart the iteration.
## 
## iterations 22
## function evaluations 55
## segments explored during Cauchy searches 27
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 0
## norm of the final projected gradient 0.00522634
## final function value -1425.92
## 
## F = -1425.92
## final  value -1425.924355 
## converged
## 
## 
## ***************************
## Estimation of CV-set 5/10
## Optimisation using starting value 1/2
## N = 5, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate     0  f=      -1452.2  |proj g|=       9.4382
## At iterate    10  f =      -1452.5  |proj g|=      0.064196
## 
## iterations 18
## function evaluations 28
## segments explored during Cauchy searches 18
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 0
## norm of the final projected gradient 0.0416277
## final function value -1452.46
## 
## F = -1452.46
## final  value -1452.455454 
## converged
## 
## Optimisation using starting value 2/2
## N = 5, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate     0  f=      -1190.7  |proj g|=           14
## At iterate    10  f =      -1452.4  |proj g|=       0.37119
## 
## iterations 19
## function evaluations 23
## segments explored during Cauchy searches 23
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 0
## norm of the final projected gradient 0.00256525
## final function value -1452.46
## 
## F = -1452.46
## final  value -1452.455454 
## converged
## 
## 
## ***************************
## Estimation of CV-set 6/10
## Optimisation using starting value 1/2
## N = 5, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate     0  f=      -1423.6  |proj g|=        10.14
## At iterate    10  f =        -1424  |proj g|=      0.020021
## 
## iterations 11
## function evaluations 29
## segments explored during Cauchy searches 12
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 0
## norm of the final projected gradient 0.0200214
## final function value -1424
## 
## F = -1424
## Warning:  more than 10 function and gradient evaluations
##    in the last line search
## final  value -1424.001300 
## converged
## 
## Optimisation using starting value 2/2
## N = 5, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate     0  f=      -1154.4  |proj g|=           14
## At iterate    10  f =        -1424  |proj g|=       0.30775
## 
## iterations 17
## function evaluations 34
## segments explored during Cauchy searches 21
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 0
## norm of the final projected gradient 0.0714014
## final function value -1424
## 
## F = -1424
## Warning:  more than 10 function and gradient evaluations
##    in the last line search
## final  value -1424.001300 
## converged
## 
## 
## ***************************
## Estimation of CV-set 7/10
## Optimisation using starting value 1/2
## N = 5, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate     0  f=      -1431.9  |proj g|=       7.9183
## At iterate    10  f =      -1434.1  |proj g|=        1.3391
## At iterate    20  f =      -1434.3  |proj g|=       0.40425
## 
## iterations 26
## function evaluations 32
## segments explored during Cauchy searches 26
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 0
## norm of the final projected gradient 0.00966739
## final function value -1434.28
## 
## F = -1434.28
## final  value -1434.275814 
## converged
## 
## Optimisation using starting value 2/2
## N = 5, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate     0  f=      -1161.7  |proj g|=           14
## At iterate    10  f =      -1434.1  |proj g|=       0.40557
## At iterate    20  f =      -1434.3  |proj g|=        0.3707
## At iterate    30  f =      -1434.3  |proj g|=      0.072941
## 
## iterations 36
## function evaluations 57
## segments explored during Cauchy searches 40
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 0
## norm of the final projected gradient 0.0147916
## final function value -1434.28
## 
## F = -1434.28
## Warning:  more than 10 function and gradient evaluations
##    in the last line search
## final  value -1434.275785 
## converged
## 
## 
## ***************************
## Estimation of CV-set 8/10
## Optimisation using starting value 1/2
## N = 5, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate     0  f=      -1423.6  |proj g|=       7.5046
## At iterate    10  f =      -1424.4  |proj g|=        2.8114
## At iterate    20  f =      -1424.6  |proj g|=       0.06844
## 
## iterations 26
## function evaluations 28
## segments explored during Cauchy searches 26
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 0
## norm of the final projected gradient 0.00782461
## final function value -1424.61
## 
## F = -1424.61
## final  value -1424.606786 
## converged
## 
## Optimisation using starting value 2/2
## N = 5, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate     0  f=      -1162.7  |proj g|=           14
## At iterate    10  f =      -1424.1  |proj g|=       0.74494
## At iterate    20  f =      -1424.6  |proj g|=      0.052911
## 
## iterations 23
## function evaluations 38
## segments explored during Cauchy searches 27
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 0
## norm of the final projected gradient 0.0144033
## final function value -1424.61
## 
## F = -1424.61
## final  value -1424.606818 
## converged
## 
## 
## ***************************
## Estimation of CV-set 9/10
## Optimisation using starting value 1/2
## N = 5, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate     0  f=      -1444.3  |proj g|=       12.803
## At iterate    10  f =        -1445  |proj g|=        0.4324
## At iterate    20  f =        -1445  |proj g|=       0.12333
## 
## iterations 21
## function evaluations 26
## segments explored during Cauchy searches 21
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 0
## norm of the final projected gradient 0.102739
## final function value -1445.01
## 
## F = -1445.01
## final  value -1445.012827 
## converged
## 
## Optimisation using starting value 2/2
## N = 5, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate     0  f=      -1174.4  |proj g|=           14
## At iterate    10  f =      -1444.9  |proj g|=        0.4452
## At iterate    20  f =        -1445  |proj g|=       0.12382
## 
## iterations 22
## function evaluations 30
## segments explored during Cauchy searches 26
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 0
## norm of the final projected gradient 0.0346693
## final function value -1445.01
## 
## F = -1445.01
## final  value -1445.012645 
## converged
## 
## 
## ***************************
## Estimation of CV-set 10/10
## Optimisation using starting value 1/2
## N = 5, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate     0  f=      -1449.9  |proj g|=        10.14
## At iterate    10  f =      -1452.2  |proj g|=        1.4276
## 
## iterations 19
## function evaluations 24
## segments explored during Cauchy searches 20
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 0
## norm of the final projected gradient 0.0118335
## final function value -1452.23
## 
## F = -1452.23
## final  value -1452.231580 
## converged
## 
## Optimisation using starting value 2/2
## N = 5, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate     0  f=        -1180  |proj g|=           14
## At iterate    10  f =      -1452.2  |proj g|=       0.75447
## At iterate    20  f =      -1452.2  |proj g|=     0.0056523
## 
## iterations 21
## function evaluations 26
## segments explored during Cauchy searches 25
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 0
## norm of the final projected gradient 0.00106247
## final function value -1452.23
## 
## F = -1452.23
## final  value -1452.231549 
## converged
## 
print(est.denver.A.cv)
## Cross-validation parameter estimation for STmodel
##   with 10 CV-groups and 2 starting points.
##   Results: 10 converged, 0 not converged.
## 
## No fixed parameters.
## 
## Estimated function values and convergence info:
##       value convergence conv  eigen.min eigen.all.min
## 1  1431.711        TRUE TRUE 0.21528208            NA
## 2  1435.925        TRUE TRUE 1.95316317            NA
## 3  1420.637        TRUE TRUE 2.37610361            NA
## 4  1425.924        TRUE TRUE 0.60717608            NA
## 5  1452.455        TRUE TRUE 1.07613640            NA
## 6  1424.001        TRUE TRUE 1.50085818            NA
## 7  1434.276        TRUE TRUE 0.06991585            NA
## 8  1424.607        TRUE TRUE 0.06886468            NA
## 9  1445.013        TRUE TRUE 0.88310030            NA
## 10 1452.232        TRUE TRUE 2.66410130            NA
par(mfrow=c(1,1), mar=c(13.5,2.5,.5,.5), las=2)
with(coef(est.denver.model.A, pars="all"),
     plotCI((1:length(par))+.3, par, uiw=1.96*sd,
             col=2, xlab="", xaxt="n", ylab=""))
boxplot(est.denver.A.cv, "all", boxwex=.4, col="grey", add=TRUE)

#' Save the results as an .rdata object
save(denver.data.A, denver.model.A, est.denver.model.A, est.denver.A.cv,
     file = here::here("Results", "Denver_ST_Model_A.rdata"))

3.1.4 Prediction using the CV model

Making predictions using the CV model. Printing out the CV summary statistics as well

pred.A.cv <- predictCV(denver.model.A, est.denver.A.cv, LTA = T)
pred.A.cv.log <- predictCV(denver.model.A, est.denver.A.cv,
                           LTA = T, transform="unbiased")

head(pred.A.cv$pred.all$EX)
##            central         d_1         d_10         d_11         d_12
## 2008-12-29      NA -0.36228648 -0.232745184 -0.216254789 -0.227981974
## 2009-01-05      NA -0.28277641 -0.161511433 -0.144077302 -0.158306997
## 2009-01-12      NA -0.12950120 -0.016765871  0.002857285 -0.023815708
## 2009-01-19      NA -0.15196738 -0.046503108 -0.027625136 -0.043730230
## 2009-01-26      NA -0.09185841  0.006385409  0.025928854  0.008920688
## 2009-02-02      NA -0.01059003  0.080960068  0.101556042  0.080178246
##                   d_13        d_14        d_15      d_16         d_17
## 2008-12-29 -0.27587718 -0.45308211 -0.41132938 0.2175431 -0.408425796
## 2009-01-05 -0.20646305 -0.37304442 -0.32350036 0.2626777 -0.318590537
## 2009-01-12 -0.07222790 -0.22837754 -0.16069099 0.3730717 -0.159188935
## 2009-01-19 -0.09238802 -0.23853533 -0.17709189 0.3300518 -0.166740871
## 2009-01-26 -0.03996632 -0.17677830 -0.10984996 0.3611382 -0.097038051
## 2009-02-02  0.03108329 -0.09725951 -0.02162496 0.4128319 -0.008492464
##                    d_18        d_19       d_2          d_20          d_21
## 2008-12-29 -0.331038019 -0.29222585 0.3901875 -0.2463607208 -0.3460726863
## 2009-01-05 -0.253982498 -0.21299537 0.4284750 -0.1730029523 -0.2677362867
## 2009-01-12 -0.103754862 -0.05749019 0.5410440 -0.0249093985 -0.1160163337
## 2009-01-19 -0.127707195 -0.08346099 0.4800610 -0.0542805594 -0.1390662727
## 2009-01-26 -0.069672425 -0.02392275 0.5039757  0.0003105882 -0.0799363876
## 2009-02-02  0.009487613  0.05771830 0.5523269  0.0768787046  0.0003005563
##                   d_22        d_23        d_25        d_26          d_28
## 2008-12-29 -0.33307875 -0.07373230 -0.25021318 -0.27685309 -0.2974751351
## 2009-01-05 -0.25255124 -0.01185879 -0.18488305 -0.20382816 -0.2294987022
## 2009-01-12 -0.09850499  0.11497190 -0.04416628 -0.05714902 -0.0861834205
## 2009-01-19 -0.11971811  0.08771213 -0.08175036 -0.08542597 -0.1212759568
## 2009-01-26 -0.05868548  0.13350760 -0.03427905 -0.03098616 -0.0714792477
## 2009-02-02  0.02333979  0.19854577  0.03600967  0.04505800  0.0009191248
##                  d_29         d_3        d_31        d_32        d_33
## 2008-12-29 -0.5944721 -0.26394191 -0.03926273 -0.44919379 -0.44143570
## 2009-01-05 -0.4996518 -0.18999990  0.01415939 -0.36925205 -0.36005273
## 2009-01-12 -0.3317455 -0.05131864  0.14181773 -0.21619031 -0.20243394
## 2009-01-19 -0.3392755 -0.06721566  0.09478918 -0.23742520 -0.22637814
## 2009-01-26 -0.2656605 -0.01081514  0.13197366 -0.17685420 -0.16494842
## 2009-02-02 -0.1722824  0.06384413  0.19244453 -0.09539324 -0.08159139
##                   d_34       d_35        d_36         d_37        d_38
## 2008-12-29 -0.27050641 -0.4494165 -0.29216521 -0.328626549 -0.28894626
## 2009-01-05 -0.19399142 -0.3648959 -0.21965675 -0.253748051 -0.21270655
## 2009-01-12 -0.04115271 -0.2107129 -0.07389406 -0.104161228 -0.06327999
## 2009-01-19 -0.06968022 -0.2232687 -0.10212755 -0.132100589 -0.08800042
## 2009-01-26 -0.01252819 -0.1582361 -0.04808845 -0.076173118 -0.03068253
## 2009-02-02  0.06694804 -0.0739274  0.02744664  0.001607336  0.04782714
##                   d_39         d_4       d_40        d_41         d_42
## 2008-12-29 -0.34337745 -0.30736245 -0.5345659 -0.36635499 -0.368793538
## 2009-01-05 -0.27222325 -0.22171269 -0.4364337 -0.29292362 -0.283314659
## 2009-01-12 -0.12715302 -0.02968994 -0.2321539 -0.14475780 -0.091459697
## 2009-01-19 -0.15748644 -0.08872076 -0.2794322 -0.17405966 -0.150651402
## 2009-01-26 -0.10472010 -0.02746285 -0.2072056 -0.11940384 -0.089543652
## 2009-02-02 -0.03011308  0.06648047 -0.1033111 -0.04277704  0.004263438
##                   d_43        d_44         d_45        d_46        d_47
## 2008-12-29 -0.35374440 -0.39385807 -0.305318718 -0.28919264 -0.27786268
## 2009-01-05 -0.27927994 -0.31137824 -0.235347298 -0.21503560 -0.20388790
## 2009-01-12 -0.13159658 -0.15541494 -0.090578943 -0.06565131 -0.05354352
## 2009-01-19 -0.15798845 -0.17478990 -0.123138425 -0.09492467 -0.08446268
## 2009-01-26 -0.10223054 -0.11204170 -0.071522991 -0.03969683 -0.02954283
## 2009-02-02 -0.02513611 -0.02846003  0.002345503  0.03762877  0.04790834
##                   d_48         d_49          d_5         d_50       d_51
## 2008-12-29 -0.28784879 -0.241872686 -0.373821107 -0.326813616 -0.5147255
## 2009-01-05 -0.21420855 -0.172601938 -0.289902488 -0.253699278 -0.4279910
## 2009-01-12 -0.06533173 -0.038507580 -0.132298341 -0.105338871 -0.2716342
## 2009-01-19 -0.09509167 -0.058802693 -0.150613781 -0.135593958 -0.2821055
## 2009-01-26 -0.04031796 -0.006506988 -0.086630858 -0.081282376 -0.2151275
## 2009-02-02  0.03659565  0.064428316 -0.001847949 -0.004788027 -0.1290538
##                   d_52        d_53        d_54       d_55         d_56
## 2008-12-29 -0.28711260 -0.12261843 -0.27544123 -0.5808377 -0.366916130
## 2009-01-05 -0.21192210 -0.05280246 -0.20218427 -0.4841291 -0.284728369
## 2009-01-12 -0.06202891  0.12367243 -0.05254475 -0.2812472 -0.128823825
## 2009-01-19 -0.08967451  0.04973375 -0.08413976 -0.3298659 -0.148768905
## 2009-01-26 -0.03347287  0.09707791 -0.02985068 -0.2588902 -0.086306954
## 2009-02-02  0.04455631  0.17839841  0.04702823 -0.1561306 -0.002903898
##                  d_57        d_58             d_6        d_60        d_61
## 2008-12-29 -0.5430766 -0.29262610 -0.323992416874 -0.30174711 -0.31409306
## 2009-01-05 -0.4525655 -0.21521568 -0.250514279586 -0.22655601 -0.23298237
## 2009-01-12 -0.2888908 -0.06440499 -0.103390117553 -0.07501728 -0.07853822
## 2009-01-19 -0.3004780 -0.08832676 -0.131240370363 -0.10479126 -0.09897610
## 2009-01-26 -0.2306498 -0.03001057 -0.076402303282 -0.04880258 -0.03740834
## 2009-02-02 -0.1407071  0.04948817  0.000003149964  0.02961824  0.04504029
##                   d_62       d_63       d_64        d_65        d_66
## 2008-12-29 -0.41729878 -0.3975971 -0.5041282 -0.39328711 -0.27193223
## 2009-01-05 -0.32935511 -0.3157525 -0.4205455 -0.30840497 -0.20124275
## 2009-01-12 -0.16802663 -0.1641974 -0.2624116 -0.15049204 -0.05662886
## 2009-01-19 -0.18225728 -0.1792727 -0.2821558 -0.16707543 -0.08739982
## 2009-01-26 -0.11470780 -0.1165916 -0.2185796 -0.10216311 -0.03504185
## 2009-02-02 -0.02677031 -0.0344163 -0.1338600 -0.01676362  0.03919469
##                   d_67        d_68         d_7         d_8 d_24 d_27 d_30 d_59
## 2008-12-29 -0.19999832 -0.35203089 -0.27310791 -0.30107370   NA   NA   NA   NA
## 2009-01-05 -0.12616073 -0.26923058 -0.19730288 -0.22400540   NA   NA   NA   NA
## 2009-01-12  0.05426315 -0.11673689 -0.04806857 -0.07714022   NA   NA   NA   NA
## 2009-01-19 -0.01588908 -0.13091238 -0.07350184 -0.09671253   NA   NA   NA   NA
## 2009-01-26  0.03498902 -0.06739144 -0.01659636 -0.03822853   NA   NA   NA   NA
## 2009-02-02  0.11951558  0.01554583  0.06162256  0.04013914   NA   NA   NA   NA
##            d_9
## 2008-12-29  NA
## 2009-01-05  NA
## 2009-01-12  NA
## 2009-01-19  NA
## 2009-01-26  NA
## 2009-02-02  NA
tail(pred.A.cv$pred.all$EX)
##            central       d_1      d_10      d_11      d_12      d_13      d_14
## 2020-11-02      NA 1.0153075 1.0365885 1.0403129 1.0482544 1.0144463 0.9451461
## 2020-11-09      NA 0.6129213 0.6313776 0.6211391 0.6343057 0.5830190 0.5403579
## 2020-11-16      NA 0.9616024 0.9988670 0.9817190 1.0003479 0.9577870 0.8917466
## 2020-11-23      NA 0.8771148 0.9064736 0.8988584 0.9086319 0.8676386 0.7912111
## 2020-11-30      NA 1.0243965 1.0570188 1.0569792 1.0663244 1.0334303 0.9434250
## 2020-12-07      NA        NA        NA        NA        NA        NA        NA
##                 d_15      d_16      d_17      d_18      d_19       d_2
## 2020-11-02 1.1085251 1.1489975 1.1710023 1.0253529 1.0865084 1.2657413
## 2020-11-09 0.7044141 0.7257865 0.7560962 0.6175884 0.6846931 0.8208796
## 2020-11-16 1.0735201 1.1230581 1.1124280 0.9765002 1.0537491 1.2079120
## 2020-11-23 0.9723052 1.0133818 1.0188022 0.8865214 0.9578493 1.1278449
## 2020-11-30 1.1125346 1.2048538 1.1681545 1.0376216 1.1021347 1.2913014
## 2020-12-07        NA        NA        NA        NA        NA        NA
##                 d_20      d_21      d_22      d_23      d_25      d_26
## 2020-11-02 1.0200958 1.0374512 1.0716942 1.0775568 0.9168185 1.0357455
## 2020-11-09 0.6135664 0.6255490 0.6636695 0.6716347 0.5135114 0.6045204
## 2020-11-16 0.9891965 0.9871666 1.0176236 1.0523952 0.8828499 0.9770651
## 2020-11-23 0.8986404 0.8977975 0.9308337 0.9539043 0.7897345 0.8918454
## 2020-11-30 1.0689908 1.0486209 1.0800445 1.1080207 0.9462193 1.0531385
## 2020-12-07        NA        NA        NA        NA        NA        NA
##                 d_28      d_29       d_3      d_31      d_32      d_33
## 2020-11-02 0.8933947 1.0406468 1.0939037 0.9419153 0.9546268 0.9500936
## 2020-11-09 0.4985135 0.6348037 0.6717188 0.5426915 0.5477225 0.5550487
## 2020-11-16 0.8561805 0.9841973 1.0426941 0.9181546 0.9064500 0.8995412
## 2020-11-23 0.7661558 0.8892249 0.9487463 0.8287602 0.8137409 0.8143944
## 2020-11-30 0.9178555 1.0339336 1.1074392 0.9826598 0.9642513 0.9577187
## 2020-12-07        NA        NA        NA        NA        NA        NA
##                 d_34      d_35      d_36      d_37      d_38      d_39
## 2020-11-02 1.0762004 1.0494578 0.9920703 0.9937540 1.0610605 0.8902132
## 2020-11-09 0.6642520 0.6338048 0.5864955 0.5850722 0.6497705 0.5031871
## 2020-11-16 1.0304952 0.9968617 0.9476558 0.9502138 1.0150518 0.8487846
## 2020-11-23 0.9387189 0.8999141 0.8576758 0.8581069 0.9195882 0.7640492
## 2020-11-30 1.0948636 1.0555747 1.0097651 1.0102803 1.0762677 0.9086520
## 2020-12-07        NA        NA        NA        NA        NA        NA
##                  d_4      d_40      d_41      d_42      d_43      d_44
## 2020-11-02 1.0395736 0.8535904 0.9384733 0.9905828 0.9672294 1.0466225
## 2020-11-09 0.6278772 0.4747605 0.5277814 0.5831185 0.5605568 0.6346899
## 2020-11-16 0.9808899 0.7888230 0.8972960 0.9403258 0.9247074 0.9905022
## 2020-11-23 0.9094250 0.6893282 0.7999222 0.8633025 0.8293803 0.9028538
## 2020-11-30 1.0531131 0.8155815 0.9586894 1.0047726 0.9844266 1.0529555
## 2020-12-07        NA        NA        NA        NA        NA        NA
##                 d_45      d_46      d_47      d_48      d_49       d_5
## 2020-11-02 0.9303928 0.9916843 1.0129575 0.9969385 1.0382311 1.0273759
## 2020-11-09 0.5281338 0.5936733 0.6169624 0.5962082 0.6248963 0.6316690
## 2020-11-16 0.8909502 0.9422248 0.9830546 0.9516216 0.9999965 0.9672184
## 2020-11-23 0.7991530 0.8602363 0.8929020 0.8655292 0.9037291 0.8789137
## 2020-11-30 0.9525905 1.0073778 1.0289339 1.0139954 1.0571023 1.0225070
## 2020-12-07        NA        NA        NA        NA        NA        NA
##                 d_50      d_51      d_52      d_53      d_54      d_55
## 2020-11-02 0.9280294 1.0091265 1.0366260 1.1157954 0.9993623 0.9486490
## 2020-11-09 0.5403765 0.6048557 0.6228188 0.6716793 0.5925683 0.5497244
## 2020-11-16 0.8832730 0.9679981 0.9849256 1.0494501 0.9637032 0.8978995
## 2020-11-23 0.8012276 0.8704512 0.8982857 0.9753801 0.8749914 0.8151933
## 2020-11-30 0.9445205 1.0101989 1.0516888 1.1248030 1.0389529 0.9508642
## 2020-12-07        NA        NA        NA        NA        NA        NA
##                 d_56      d_57      d_58       d_6      d_60      d_61
## 2020-11-02 1.0291476 1.0143207 1.0992731 0.9797674 0.9921893 1.0747313
## 2020-11-09 0.6269554 0.6146639 0.6820153 0.5707252 0.5870810 0.6667059
## 2020-11-16 0.9685909 0.9605498 1.0482364 0.9351738 0.9428642 1.0134782
## 2020-11-23 0.8847239 0.8682101 0.9582192 0.8446043 0.8522978 0.9295177
## 2020-11-30 1.0316212 1.0118100 1.1078981 0.9969283 1.0040087 1.0795107
## 2020-12-07        NA        NA        NA        NA        NA        NA
##                 d_62      d_63      d_64      d_65      d_66      d_67
## 2020-11-02 1.0705972 1.0484585 0.9326397 1.0942474 0.9797545 0.9938950
## 2020-11-09 0.6681986 0.6388058 0.5362505 0.6543527 0.5715386 0.5896571
## 2020-11-16 1.0063436 0.9954318 0.8829836 1.0208160 0.9361592 0.9660635
## 2020-11-23 0.9218063 0.9041185 0.7933088 0.9373491 0.8440889 0.8821866
## 2020-11-30 1.0666864 1.0544009 0.9377339 1.1003707 1.0005123 1.0210878
## 2020-12-07        NA        NA        NA        NA        NA        NA
##                 d_68       d_7       d_8 d_24 d_27 d_30 d_59 d_9
## 2020-11-02 1.1566223 1.0348270 1.0718938   NA   NA   NA   NA  NA
## 2020-11-09 0.7310982 0.6245033 0.6607406   NA   NA   NA   NA  NA
## 2020-11-16 1.0958777 0.9774879 1.0209648   NA   NA   NA   NA  NA
## 2020-11-23 1.0038295 0.8910882 0.9309178   NA   NA   NA   NA  NA
## 2020-11-30 1.1551023 1.0442159 1.0832408   NA   NA   NA   NA  NA
## 2020-12-07        NA        NA        NA   NA   NA   NA   NA  NA
head(pred.A.cv.log$pred.all$EX)
##            central       d_1      d_10      d_11      d_12      d_13      d_14
## 2008-12-29      NA 0.7142497 0.8123136 0.8269625 0.8179457 0.7796934 0.6530775
## 2009-01-05      NA 0.7732058 0.8722277 0.8886006 0.8766751 0.8354583 0.7072605
## 2009-01-12      NA 0.9011275 1.0080113 1.0289849 1.0025813 0.9552003 0.8171081
## 2009-01-19      NA 0.8809779 0.9784260 0.9978760 0.9825665 0.9359015 0.8086473
## 2009-01-26      NA 0.9354427 1.0315220 1.0525857 1.0354722 0.9860685 0.8599840
## 2009-02-02      NA 1.0145411 1.1113505 1.1351168 1.1117679 1.0585039 0.9310091
##                 d_15     d_16      d_17      d_18      d_19      d_2      d_20
## 2008-12-29 0.6793596 1.277067 0.6818200 0.7365967 0.7656016 1.518451 0.8024370
## 2009-01-05 0.7417020 1.335580 0.7457778 0.7954396 0.8286148 1.577126 0.8632654
## 2009-01-12 0.8728193 1.491027 0.8745214 0.9242156 0.9679108 1.764458 1.0008064
## 2009-01-19 0.8586003 1.427886 0.8678293 0.9022054 0.9430001 1.659605 0.9716286
## 2009-01-26 0.9183014 1.472667 0.9303781 0.9559953 1.0007652 1.699380 1.0259628
## 2009-02-02 1.0029841 1.550543 1.0164302 1.0346460 1.0858237 1.783240 1.1074480
##                 d_21      d_22      d_23      d_25      d_26      d_28
## 2008-12-29 0.7252808 0.7366993 0.9543645 0.7981263 0.7793073 0.7612828
## 2009-01-05 0.7843267 0.7981790 1.0149399 0.8519815 0.8380318 0.8148053
## 2009-01-12 0.9127702 0.9308048 1.1518457 0.9806869 0.9701062 0.9403349
## 2009-01-19 0.8919249 0.9110118 1.1205897 0.9444900 0.9427942 0.9078866
## 2009-01-26 0.9462141 0.9681217 1.1728591 0.9903878 0.9953129 0.9542220
## 2009-02-02 1.0252295 1.0506875 1.2514702 1.0624890 1.0737561 1.0258522
##                 d_29       d_3      d_31      d_32      d_33      d_34
## 2008-12-29 0.5657541 0.7890550 0.9865918 0.6545086 0.6594803 0.7824119
## 2009-01-05 0.6219814 0.8493264 1.0405254 0.7088368 0.7152971 0.8445123
## 2009-01-12 0.7356504 0.9753831 1.1820029 0.8259299 0.8373115 0.9838539
## 2009-01-19 0.7300936 0.9597593 1.1275350 0.8084543 0.8174169 0.9560853
## 2009-01-26 0.7858334 1.0152366 1.1701087 0.8588302 0.8691315 1.0122337
## 2009-02-02 0.8627188 1.0937557 1.2429293 0.9316282 0.9446216 1.0958920
##                 d_35      d_36      d_37      d_38      d_39       d_4
## 2008-12-29 0.6544368 0.7657941 0.7390663 0.7682631 0.7278840 0.7547770
## 2009-01-05 0.7120317 0.8232176 0.7963009 0.8289591 0.7814088 0.8220648
## 2009-01-12 0.8306037 0.9522296 0.9245522 0.9623904 0.9032459 0.9958726
## 2009-01-19 0.8201335 0.9255812 0.8988836 0.9387498 0.8761291 0.9386090
## 2009-01-26 0.8751480 0.9768538 0.9504191 0.9940057 0.9234883 0.9977483
## 2009-02-02 0.9520495 1.0533951 1.0271489 1.0750841 0.9949263 1.0958884
##                 d_40      d_41      d_42      d_43      d_44      d_45
## 2008-12-29 0.6013752 0.7117019 0.7098057 0.7200597 0.6932568 0.7564946
## 2009-01-05 0.6632144 0.7657085 0.7729522 0.7755694 0.7525790 0.8110891
## 2009-01-12 0.8133450 0.8877702 0.9362192 0.8988387 0.8793119 0.9371954
## 2009-01-19 0.7756392 0.8619476 0.8822437 0.8752950 0.8621972 0.9069758
## 2009-01-26 0.8336039 0.9102072 0.9376908 0.9253712 0.9178203 0.9548490
## 2009-02-02 0.9247552 0.9825564 1.0297833 0.9994358 0.9976476 1.0279073
##                 d_46      d_47      d_48      d_49       d_5      d_50
## 2008-12-29 0.7676143 0.7766774 0.7686466 0.8066624 0.7060585 0.7392724
## 2009-01-05 0.8266755 0.8361958 0.8273595 0.8642322 0.7677155 0.7953232
## 2009-01-12 0.9598415 0.9717383 0.9601483 0.9879592 0.8986104 0.9224938
## 2009-01-19 0.9321286 0.9420559 0.9319729 0.9678678 0.8821712 0.8949801
## 2009-01-26 0.9850366 0.9951566 0.9844249 1.0196199 0.9403456 0.9449133
## 2009-02-02 1.0642106 1.0752239 1.0631117 1.0943948 1.0234492 1.0200142
##                 d_51      d_52      d_53      d_54      d_55      d_56
## 2008-12-29 0.6130619 0.7703936 0.9079292 0.7785604 0.5741824 0.7109506
## 2009-01-05 0.6684939 0.8303133 0.9733361 0.8376216 0.6323246 0.7716980
## 2009-01-12 0.7815128 0.9643380 1.1609358 0.9727094 0.7743795 0.9017381
## 2009-01-19 0.7732716 0.9378403 1.0779900 0.9423602 0.7374910 0.8838002
## 2009-01-26 0.8267494 0.9918811 1.1300779 0.9948502 0.7916139 0.9406502
## 2009-02-02 0.9009868 1.0722249 1.2256647 1.0742780 0.8771777 1.0223690
##                 d_57      d_58       d_6      d_60      d_61      d_62
## 2008-12-29 0.5955915 0.7650992 0.7434237 0.7583467 0.7488499 0.6771953
## 2009-01-05 0.6519686 0.8266209 0.7998065 0.8174541 0.8120643 0.7391708
## 2009-01-12 0.7678618 0.9611164 0.9262688 0.9510939 0.9476283 0.8682919
## 2009-01-19 0.7589761 0.9383486 0.9005751 0.9230986 0.9284088 0.8557829
## 2009-01-26 0.8138332 0.9946538 0.9511207 0.9761735 0.9873227 0.9153765
## 2009-02-02 0.8903940 1.0769189 1.0264518 1.0557367 1.0721395 0.9993348
##                 d_63      d_64      d_65      d_66      d_67      d_68
## 2008-12-29 0.6892434 0.6201032 0.6921422 0.7817906 0.8403231 0.7213761
## 2009-01-05 0.7478973 0.6739660 0.7533067 0.8388895 0.9044899 0.7835130
## 2009-01-12 0.8701524 0.7892311 0.8820142 0.9692466 1.0830890 0.9124459
## 2009-01-19 0.8570216 0.7736336 0.8673772 0.9397370 1.0095204 0.8994860
## 2009-01-26 0.9123626 0.8242687 0.9254336 0.9901302 1.0620464 0.9583736
## 2009-02-02 0.9904190 0.8970170 1.0078387 1.0663282 1.1555778 1.0411594
##                  d_7       d_8 d_24 d_27 d_30 d_59 d_9
## 2008-12-29 0.7801793 0.7590881   NA   NA   NA   NA  NA
## 2009-01-05 0.8415614 0.8197610   NA   NA   NA   NA  NA
## 2009-01-12 0.9769465 0.9493006   NA   NA   NA   NA  NA
## 2009-01-19 0.9523631 0.9307804   NA   NA   NA   NA  NA
## 2009-01-26 1.0080862 0.9867341   NA   NA   NA   NA  NA
## 2009-02-02 1.0900662 1.0670824   NA   NA   NA   NA  NA
tail(pred.A.cv.log$pred.all$EX)
##            central      d_1     d_10     d_11     d_12     d_13     d_14
## 2020-11-02      NA 2.781967 2.834168 2.850092 2.871261 2.773636 2.590173
## 2020-11-09      NA 1.860271 1.889887 1.874058 1.897848 1.801553 1.727810
## 2020-11-16      NA 2.636266 2.729150 2.687548 2.736549 2.620466 2.455126
## 2020-11-23      NA 2.422608 2.488260 2.473722 2.496600 2.394446 2.220186
## 2020-11-30      NA 2.806957 2.892492 2.897385 2.922914 2.826112 2.585098
## 2020-12-07      NA       NA       NA       NA       NA       NA       NA
##                d_15     d_16     d_17     d_18     d_19      d_2     d_20
## 2020-11-02 3.042630 3.168553 3.244899 2.804743 2.977606 3.565357 2.787609
## 2020-11-09 2.031147 2.075041 2.142842 1.865439 1.992255 2.284875 1.856295
## 2020-11-16 2.937909 3.086961 3.060044 2.670773 2.881450 3.364473 2.702442
## 2020-11-23 2.655084 2.766155 2.786469 2.440876 2.617895 3.105413 2.468360
## 2020-11-30 3.054758 3.349759 3.235245 2.838943 3.024171 3.656686 2.926678
## 2020-12-07       NA       NA       NA       NA       NA       NA       NA
##                d_21     d_22     d_23     d_25     d_26     d_28     d_29
## 2020-11-02 2.838464 2.944943 2.950946 2.514945 2.836488 2.459586 2.849619
## 2020-11-09 1.880129 1.958100 1.966235 1.680234 1.842735 1.657154 1.898990
## 2020-11-16 2.699164 2.789486 2.877194 2.430905 2.674396 2.369693 2.693118
## 2020-11-23 2.468375 2.557450 2.607189 2.214754 2.455786 2.165668 2.449085
## 2020-11-30 2.870175 2.968854 3.041489 2.589904 2.885488 2.520415 2.830385
## 2020-12-07       NA       NA       NA       NA       NA       NA       NA
##                 d_3     d_31     d_32     d_33     d_34     d_35     d_36
## 2020-11-02 3.003200 2.578411 2.612859 2.606174 2.949157 2.871160 2.712563
## 2020-11-09 1.968774 1.729614 1.739312 1.755585 1.953326 1.894619 1.808085
## 2020-11-16 2.852857 2.517638 2.489736 2.477540 2.817212 2.723829 2.594486
## 2020-11-23 2.596908 2.302268 2.269219 2.275265 2.570108 2.472086 2.371153
## 2020-11-30 3.043392 2.685236 2.637735 2.625844 3.004389 2.888393 2.760578
## 2020-12-07       NA       NA       NA       NA       NA       NA       NA
##                d_37     d_38     d_39      d_4     d_40     d_41     d_42
## 2020-11-02 2.717994 2.903920 2.455849 2.851992 2.365285 2.570651 2.713591
## 2020-11-09 1.806047 1.924605 1.667618 1.889399 1.619319 1.704710 1.805334
## 2020-11-16 2.601857 2.773088 2.355974 2.689141 2.216697 2.466630 2.580294
## 2020-11-23 2.372806 2.520521 2.164496 2.503570 2.006685 2.237667 2.388918
## 2020-11-30 2.762705 2.947977 2.501186 2.890340 2.276651 2.622596 2.751865
## 2020-12-07       NA       NA       NA       NA       NA       NA       NA
##                d_43     d_44     d_45     d_46     d_47     d_48     d_49
## 2020-11-02 2.644784 2.871290 2.552176 2.716652 2.767308 2.728546 2.838609
## 2020-11-09 1.760972 1.901682 1.706792 1.824634 1.862357 1.827645 1.877416
## 2020-11-16 2.534449 2.714152 2.453156 2.585510 2.685604 2.607608 2.731720
## 2020-11-23 2.303931 2.486247 2.237889 2.381970 2.454025 2.392490 2.480876
## 2020-11-30 2.690259 2.888769 2.608915 2.759540 2.811569 2.775402 2.891986
## 2020-12-07       NA       NA       NA       NA       NA       NA       NA
##                 d_5     d_50     d_51     d_52     d_53     d_54     d_55
## 2020-11-02 2.816501 2.550457 2.756640 2.838729 3.068558 2.730138 2.603335
## 2020-11-09 1.895985 1.730845 1.839873 1.876630 1.968020 1.817608 1.746836
## 2020-11-16 2.651826 2.438780 2.645348 2.695350 2.871252 2.634324 2.474236
## 2020-11-23 2.427622 2.246664 2.399420 2.471544 2.666157 2.410638 2.277745
## 2020-11-30 2.802412 2.592789 2.759231 2.881209 3.095746 2.840086 2.608630
## 2020-12-07       NA       NA       NA       NA       NA       NA       NA
##                d_56     d_57     d_58      d_6     d_60     d_61     d_62
## 2020-11-02 2.821899 2.777065 3.017327 2.683165 2.713423 2.950928 2.945523
## 2020-11-09 1.887339 1.862124 1.987928 1.782228 1.809528 1.962216 1.969537
## 2020-11-16 2.655848 2.631589 2.867090 2.565725 2.582660 2.775499 2.761772
## 2020-11-23 2.442118 2.399440 2.620244 2.343427 2.358982 2.551949 2.537752
## 2020-11-30 2.828475 2.769939 3.043280 2.728884 2.745389 2.964892 2.933258
## 2020-12-07       NA       NA       NA       NA       NA       NA       NA
##                d_63     d_64     d_65     d_66     d_67     d_68      d_7
## 2020-11-02 2.871315 2.561443 3.004188 2.679337 2.715772 3.195663 2.833098
## 2020-11-09 1.906125 1.723075 1.934912 1.781227 1.812622 2.088038 1.879540
## 2020-11-16 2.722806 2.437048 2.791237 2.564807 2.640928 3.007077 2.675123
## 2020-11-23 2.485121 2.227918 2.567636 2.339134 2.428353 2.742560 2.453666
## 2020-11-30 2.888048 2.573989 3.022190 2.735129 2.790114 3.190385 2.859654
## 2020-12-07       NA       NA       NA       NA       NA       NA       NA
##                 d_8 d_24 d_27 d_30 d_59 d_9
## 2020-11-02 2.938774   NA   NA   NA   NA  NA
## 2020-11-09 1.947983   NA   NA   NA   NA  NA
## 2020-11-16 2.792628   NA   NA   NA   NA  NA
## 2020-11-23 2.552077   NA   NA   NA   NA  NA
## 2020-11-30 2.971919   NA   NA   NA   NA  NA
## 2020-12-07       NA   NA   NA   NA   NA  NA
summary(pred.A.cv)
## Cross-validation predictions for STmodel with 10 CV-groups.
##   Predictions for 790 observations.
##   Temporal averages for 63 locations.
## 
## RMSE:
##              EX.mu EX.mu.beta         EX
## obs     0.16756831 0.16756831 0.10939771
## average 0.09377047 0.09377047 0.05212872
## 
## R2:
##             EX.mu EX.mu.beta        EX
## obs     0.5831563  0.5831563 0.8223333
## average 0.6093873  0.6093873 0.8792830
## 
## Coverage of 95% prediction intervals:
##                EX
## obs     0.9443038
## average 0.8571429
summary(pred.A.cv.log)
## Cross-validation predictions for STmodel with 10 CV-groups.
##   Predictions for 790 observations.
##   Temporal averages for 63 locations.
## 
## RMSE:
##             EX.mu EX.mu.beta         EX    EX.pred
## obs     0.2182824  0.2182824 0.15097844 0.15098707
## average 0.1045726  0.1045726 0.06003835 0.06050165
## 
## R2:
##             EX.mu EX.mu.beta        EX   EX.pred
## obs     0.6108816  0.6108816 0.8138452 0.8138239
## average 0.6750348  0.6750348 0.8928829 0.8912233
## 
## Coverage of 95% prediction intervals:
##           EX.pred
## obs     0.9443038
## average 0.9047619
par(mfrow=c(1,2), mar=c(3.3,3.3,1.5,1), mgp=c(2,1,0))
plot(pred.A.cv, "obs", ID="all", pch=c(19,NA), cex=.25, lty=c(NA,2),
     col=c("ID", "black", "grey"),
     ylim=c(-1,2),
     xlab="Observations", ylab="Predictions",
     main="Cross-validation BC (log ug/m3)")
with(pred.A.cv.log$pred.LTA, plotCI(obs, EX.pred, uiw=1.96*sqrt(VX.pred),
                                    xlab="Observations", ylab="Predictions",
                                    main="Temporal average BC (ug/m3)"))
abline(0, 1, col="grey")

jpeg(filename = here::here("Figs", "ST_CV_Obs_vs_Pred_BC_ModA.jpeg"),
     width = 8, height = 4, units = "in", res = 500)
par(mfrow=c(1,2), mar=c(3.3,3.3,1.5,1), mgp=c(2,1,0))
plot(pred.A.cv, "obs", ID="all", pch=c(19,NA), cex=.25, lty=c(NA,2),
     col=c("ID", "black", "grey"),
     ylim=c(-1,2),
     xlab="Observations", ylab="Predictions",
     main="Cross-validation BC (log ug/m3)")
with(pred.A.cv.log$pred.LTA, plotCI(obs, EX.pred, uiw=1.96*sqrt(VX.pred),
                                      xlab="Observations", ylab="Predictions",
                                      main="Temporal average BC (ug/m3)"))
abline(0, 1, col="grey")
dev.off()
## quartz_off_screen 
##                 2

3.1.5 Predictions for 2009-2020

What do the long-term predictions look like for this model? Predicting at the distributed (residential + community) sites

x.A <- coef(est.denver.model.A, pars = "cov")$par
x.A
## [1] -7.797037 -7.526946 12.449157 -3.008551 -4.860112
E.A <- predict(denver.model.A, est.denver.model.A, LTA = T,
               transform="unbiased", pred.var=FALSE)
print(E.A)
## Prediction for STmodel.
## 
## Regression parameters:
##   0 Spatio-temporal covariate(s).
##   28 beta-fields regression parameters in x$pars.
## 
## Regression parameters are assumed to be known and
##  prediction variances do NOT include
##  uncertainties in regression parameters.
## 
## Prediction of beta-fields, (x$beta):
## List of 2
##  $ mu: num [1:69, 1:2] 0.2557 0.028 0.0441 0.0712 0.042 ...
##   ..- attr(*, "dimnames")=List of 2
##  $ EX: num [1:69, 1:2] 0.2646 0.0409 0.0405 0.062 0.0302 ...
##   ..- attr(*, "dimnames")=List of 2
## 
## Predictions for log-Gaussian field of type: unbiased 
## 
## Predictions for 622 times at 69 locations.
## List of 5
##  $ EX.mu     : num [1:622, 1:69] 1.21 1.27 1.43 1.36 1.41 ...
##   ..- attr(*, "dimnames")=List of 2
##  $ EX.mu.beta: num [1:622, 1:69] 1.26 1.32 1.48 1.41 1.45 ...
##   ..- attr(*, "dimnames")=List of 2
##  $ EX        : num [1:622, 1:69] 1.29 1.35 1.52 1.44 1.48 ...
##   ..- attr(*, "dimnames")=List of 2
##  $ EX.pred   : num [1:622, 1:69] 1.3 1.36 1.53 1.45 1.49 ...
##   ..- attr(*, "dimnames")=List of 2
##  $ log.EX    : num [1:622, 1:69] 0.231 0.276 0.394 0.341 0.37 ...
##   ..- attr(*, "dimnames")=List of 2
## 
## Variances have been computed.
## List of 2
##  $ log.VX     : num [1:622, 1:69] 0.0499 0.0498 0.0497 0.0497 0.0496 ...
##   ..- attr(*, "dimnames")=List of 2
##  $ log.VX.pred: num [1:622, 1:69] 0.0576 0.0575 0.0575 0.0574 0.0574 ...
##   ..- attr(*, "dimnames")=List of 2
## 
## Mean squared prediciton errors NOT been computed.
## 
## 69 temporal averages have been compute.
## List of 1
##  $ LTA:'data.frame': 69 obs. of  4 variables:
week_preds <- as.data.frame(E.A$EX) %>%
  mutate(week = as.Date(rownames(E.A$EX)),
         year = year(as.Date(rownames(E.A$EX))))

week_sites <- colnames(week_preds)[which(str_detect(colnames(week_preds), "d_"))]

box_preds <- week_preds %>%
  select(week, all_of(week_sites)) %>%
  #filter(week %in% week_weeks) %>%
  mutate(month = month(week),
         year = year(week)) %>%
  filter(year %in% c(2009:2020))

box_data <- box_preds %>%
  pivot_longer(-c(week, month, year), names_to = "location", values_to = "pred")
summary(box_data)
##       week                month             year        location        
##  Min.   :2009-01-05   Min.   : 1.000   Min.   :2009   Length:42228      
##  1st Qu.:2012-01-09   1st Qu.: 4.000   1st Qu.:2012   Class :character  
##  Median :2014-12-29   Median : 7.000   Median :2014   Mode  :character  
##  Mean   :2014-12-28   Mean   : 6.494   Mean   :2014                     
##  3rd Qu.:2017-12-18   3rd Qu.: 9.000   3rd Qu.:2017                     
##  Max.   :2020-12-07   Max.   :12.000   Max.   :2020                     
##                                                                         
##       pred       
##  Min.   :0.3202  
##  1st Qu.:1.1202  
##  Median :1.3738  
##  Mean   :1.4175  
##  3rd Qu.:1.6484  
##  Max.   :3.8293  
##  NA's   :861
#' Box plot of weekly BC predicted at all distributed sites grouped by month
box_summary <- box_data %>%
  group_by(month) %>%
  summarize(median = median(pred, na.rm=T)) %>%
  arrange(desc(median))
box_summary
pred_box_plot <- ggplot(box_data) +
  geom_boxplot(aes(x = as.factor(month), y = pred), #, color = as.factor(month)),
               show.legend = F) +
  #scale_color_viridis(name = "Month", discrete = T) +
  facet_wrap(. ~ year, scales = "free") +
  xlab("") + ylab("Distributed site BC concentration (\u03BCg/m\u00B3): Predicted") +
  # theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust = 1)) +
  scale_x_discrete(labels = str_sub(month.name, 1, 1))
pred_box_plot

3.2 Model B: 2 temporal trend (basis) functions

3.2.1 Create the model object

For this version of the model, use iid for cov.beta (beta, beta1, beta2) and exp for cov.nu (error). Here we can specify different LUR formluae. The length of the LUR list should be number of basis functions + 1.

denver.data.B <- denver.data2.2
names(denver.data.B$covars)
##  [1] "ID"              "lon"             "lat"             "tree_cover_100" 
##  [5] "impervious_2500" "low_int_100"     "med_int_50"      "med_int_2500"   
##  [9] "high_int_100"    "pop_den_50"      "dist_m_cafos"    "dist_m_og_wells"
## [13] "dist_m_wwtp"     "aadt_100"        "aadt_250"        "aadt_1000"      
## [17] "x"               "y"               "type"
LUR.B <- list(covar_fun, covar_fun, covar_fun)

cov.beta.B <-  list(covf = c("iid", "iid", "iid"), nugget = T)
cov.nu.B <- list(covf = "exp", nugget = T, random.effect = FALSE)
locations.B <- list(coords = c("x","y"), long.lat = c("lon","lat"), others= "type")

denver.model.B <- createSTmodel(denver.data.B, LUR = LUR.B,
                                ST = c("bc_st_no2", "bc_st_smk"),
                                cov.beta = cov.beta.B, cov.nu = cov.nu.B,
                                locations = locations.B)
## No trend $trend.fnc object detected, STdata probably from old version of the package.
## $trend.fnc has been added based on spline fit to elements in STmodel$trend.
denver.model.B
## STmodel-object with:
##  No. locations: 69 (observed: 64)
##  No. time points: 622 (observed: 231)
##  No. obs: 1021
## 
## Trend with 2 basis function(s):
## [1] "V1" "V2"
## with dates:
##  2008-12-29 to 2020-12-07
## 
## Models for the beta-fields are:
## $const
## ~~tree_cover_100 + impervious_2500 + low_int_100 + med_int_50 + 
##     med_int_2500 + high_int_100 + pop_den_50 + dist_m_cafos + 
##     dist_m_og_wells + dist_m_wwtp + aadt_100 + aadt_250 + aadt_1000
## 
## $V1
## ~~tree_cover_100 + impervious_2500 + low_int_100 + med_int_50 + 
##     med_int_2500 + high_int_100 + pop_den_50 + dist_m_cafos + 
##     dist_m_og_wells + dist_m_wwtp + aadt_100 + aadt_250 + aadt_1000
## 
## $V2
## ~~tree_cover_100 + impervious_2500 + low_int_100 + med_int_50 + 
##     med_int_2500 + high_int_100 + pop_den_50 + dist_m_cafos + 
##     dist_m_og_wells + dist_m_wwtp + aadt_100 + aadt_250 + aadt_1000
## 
## 2 spatio-temporal covariate(s):
## [1] "bc_st_no2" "bc_st_smk"
## 
## Covariance model for the beta-field(s):
##  Covariance type(s): iid, iid, iid 
##  Nugget: Yes, Yes, Yes 
## Covariance model for the nu-field(s):
##  Covariance type: exp 
##  Nugget: ~1 
##  Random effect: No 
## All sites:
## central    dist 
##       1      68 
## Observed:
## central    dist 
##       1      63 
## 
## For central:
##   Number of obs: 231
##   Dates: 2016-02-08 to 2020-11-30
## For dist:
##   Number of obs: 790
##   Dates: 2018-05-07 to 2020-04-06

3.2.2 Estimate model parameters

  • nugget values can range from -1 to -6
  • range values can range from 0 to 4
  • try starting values where the sill and nugget are the same (e.g., both 0) and where they are very different (e.g., 0 and -5)
  • make sure the starting values are pretty different
names <- loglikeSTnames(denver.model.B, all=FALSE)
names
## [1] "log.nugget.const.iid"          "log.nugget.V1.iid"            
## [3] "log.nugget.V2.iid"             "nu.log.range.exp"             
## [5] "nu.log.sill.exp"               "nu.log.nugget.(Intercept).exp"
# x.init.B <- cbind(c(0, 0, 0, 0, 0, 0),
#                   c(-1, -1, -1, 0, -1, -1),
#                   c(-1, -1, -1, 0, -5, -1),
#                   c(-5, -5, -5, 0, -1, -5),
#                   c(-5, -5, -5, 0, -5, -5),
#                   c(-1, -1, -1, 2, -1, -1),
#                   c(-1, -1, -1, 2, -5, -1),
#                   c(-5, -5, -5, 2, -1, -5),
#                   c(-5, -5, -5, 2, -5, -5),
#                   c(-1, -1, -1, 4, -1, -1),
#                   c(-1, -1, -1, 4, -5, -1),
#                   c(-5, -5, -5, 4, -1, -5),
#                   c(-5, -5, -5, 4, -5, -5))

x.init.B <- cbind(c(0, 0, 0, 0, 0, 0),
                  c(-10, -5, -5, 4, -3, -5),
                  c(-10, -5, -5, 4, -5, -5),
                  c(-10, -5, -5, 6, -3, -5),
                  c(-10, -5, -5, 6, -5, -5),
                  c(-10, -5, -5, 8, -3, -5),
                  c(-10, -5, -5, 8, -5, -5),
                  c(-10, -5, -5, 10, -3, -5),
                  c(-10, -5, -5, 10, -5, -5),
                  c(-12, -5, -5, 8, -5, -5),
                  c(-12, -5, -5, 10, -5, -5),
                  c(-12, -5, -5, 12, -5, -5))

rownames(x.init.B) <- loglikeSTnames(denver.model.B, all=FALSE)
x.init.B
##                               [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9]
## log.nugget.const.iid             0  -10  -10  -10  -10  -10  -10  -10  -10
## log.nugget.V1.iid                0   -5   -5   -5   -5   -5   -5   -5   -5
## log.nugget.V2.iid                0   -5   -5   -5   -5   -5   -5   -5   -5
## nu.log.range.exp                 0    4    4    6    6    8    8   10   10
## nu.log.sill.exp                  0   -3   -5   -3   -5   -3   -5   -3   -5
## nu.log.nugget.(Intercept).exp    0   -5   -5   -5   -5   -5   -5   -5   -5
##                               [,10] [,11] [,12]
## log.nugget.const.iid            -12   -12   -12
## log.nugget.V1.iid                -5    -5    -5
## log.nugget.V2.iid                -5    -5    -5
## nu.log.range.exp                  8    10    12
## nu.log.sill.exp                  -5    -5    -5
## nu.log.nugget.(Intercept).exp    -5    -5    -5
#' Difference from tutorial: use Josh Keller's version of the function
source(here::here("Code", "functions_model.R"))
est.denver.model.B <- estimate.STmodel(denver.model.B, x.init.B)
## Optimisation using starting value 1/12
## N = 6, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate     0  f=       507.58  |proj g|=           15
## At iterate    10  f =      -1393.9  |proj g|=        2.5132
## At iterate    20  f =      -1394.2  |proj g|=      0.013016
## 
## iterations 21
## function evaluations 35
## segments explored during Cauchy searches 26
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 1
## norm of the final projected gradient 0.0130169
## final function value -1394.21
## 
## F = -1394.21
## final  value -1394.206501 
## converged
## Optimisation using starting value 2/12
## N = 6, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate     0  f=      -1197.7  |proj g|=           12
## At iterate    10  f =      -1393.8  |proj g|=        3.6232
## At iterate    20  f =      -1394.1  |proj g|=        2.0253
## ys=-8.565e+00  -gs= 2.159e+00, BFGS update SKIPPED
## At iterate    30  f =      -1608.5  |proj g|=        20.131
## 
## iterations 38
## function evaluations 58
## segments explored during Cauchy searches 41
## BFGS updates skipped 1
## active bounds at final generalized Cauchy point 0
## norm of the final projected gradient 0.168749
## final function value -1622.14
## 
## F = -1622.14
## final  value -1622.137478 
## converged
## Optimisation using starting value 3/12
## N = 6, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate     0  f=      -1237.6  |proj g|=           20
## At iterate    10  f =      -1392.3  |proj g|=        3.6831
## At iterate    20  f =      -1394.1  |proj g|=       0.28437
## ys=-5.332e-02  -gs= 1.204e+00, BFGS update SKIPPED
## At iterate    30  f =      -1606.3  |proj g|=        18.399
## At iterate    40  f =      -1622.1  |proj g|=       0.40807
## 
## iterations 45
## function evaluations 68
## segments explored during Cauchy searches 49
## BFGS updates skipped 1
## active bounds at final generalized Cauchy point 0
## norm of the final projected gradient 0.0219034
## final function value -1622.13
## 
## F = -1622.13
## final  value -1622.130756 
## converged
## Optimisation using starting value 4/12
## N = 6, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate     0  f=      -1202.8  |proj g|=           12
## At iterate    10  f =      -1631.2  |proj g|=        7.6237
## At iterate    20  f =      -1635.5  |proj g|=       0.79559
## At iterate    30  f =      -1635.5  |proj g|=       0.11189
## At iterate    40  f =      -1635.5  |proj g|=      0.090055
## 
## iterations 43
## function evaluations 62
## segments explored during Cauchy searches 47
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 0
## norm of the final projected gradient 0.00445772
## final function value -1635.49
## 
## F = -1635.49
## final  value -1635.494460 
## converged
## Optimisation using starting value 5/12
## N = 6, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate     0  f=      -1245.1  |proj g|=           20
## At iterate    10  f =      -1613.5  |proj g|=        3.2561
## At iterate    20  f =      -1624.6  |proj g|=        7.3182
## At iterate    30  f =      -1626.9  |proj g|=        4.5415
## At iterate    40  f =      -1630.8  |proj g|=        1.2802
## At iterate    50  f =      -1631.1  |proj g|=         1.002
## ys=-4.238e-01  -gs= 3.760e-01, BFGS update SKIPPED
## At iterate    60  f =      -1635.3  |proj g|=        1.4564
## At iterate    70  f =      -1635.5  |proj g|=       0.35894
## At iterate    80  f =      -1635.5  |proj g|=        0.2711
## 
## iterations 88
## function evaluations 136
## segments explored during Cauchy searches 93
## BFGS updates skipped 1
## active bounds at final generalized Cauchy point 0
## norm of the final projected gradient 0.0246402
## final function value -1635.49
## 
## F = -1635.49
## Warning:  more than 10 function and gradient evaluations
##    in the last line search
## final  value -1635.494327 
## converged
## Optimisation using starting value 6/12
## N = 6, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate     0  f=      -1262.8  |proj g|=           12
## At iterate    10  f =      -1635.3  |proj g|=        3.9677
## At iterate    20  f =      -1635.5  |proj g|=      0.072351
## Bad direction in the line search;
##    refresh the lbfgs memory and restart the iteration.
## Line search cannot locate an adequate point after 20 function
## and gradient evaluations
## final  value -1635.492358 
## stopped after 28 iterations
## Optimisation using starting value 7/12
## N = 6, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate     0  f=      -1297.5  |proj g|=           20
## At iterate    10  f =      -1634.5  |proj g|=        2.4806
## At iterate    20  f =      -1635.5  |proj g|=      0.011806
## At iterate    30  f =      -1635.5  |proj g|=      0.063142
## At iterate    40  f =      -1635.5  |proj g|=     0.0048559
## 
## iterations 41
## function evaluations 69
## segments explored during Cauchy searches 46
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 0
## norm of the final projected gradient 0.00485598
## final function value -1635.49
## 
## F = -1635.49
## Warning:  more than 10 function and gradient evaluations
##    in the last line search
## final  value -1635.494445 
## converged
## Optimisation using starting value 8/12
## N = 6, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate     0  f=      -1459.9  |proj g|=           12
## At iterate    10  f =      -1635.4  |proj g|=        2.6314
## At iterate    20  f =      -1635.5  |proj g|=      0.029661
## 
## iterations 22
## function evaluations 26
## segments explored during Cauchy searches 27
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 0
## norm of the final projected gradient 0.00905988
## final function value -1635.49
## 
## F = -1635.49
## final  value -1635.492637 
## converged
## Optimisation using starting value 9/12
## N = 6, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate     0  f=      -1370.4  |proj g|=           20
## At iterate    10  f =      -1635.5  |proj g|=       0.85406
## Bad direction in the line search;
##    refresh the lbfgs memory and restart the iteration.
## Line search cannot locate an adequate point after 20 function
## and gradient evaluations
## final  value -1635.492327 
## stopped after 14 iterations
## Optimisation using starting value 10/12
## N = 6, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate     0  f=      -1297.3  |proj g|=           20
## At iterate    10  f =      -1634.9  |proj g|=         2.363
## 
## iterations 19
## function evaluations 42
## segments explored during Cauchy searches 24
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 0
## norm of the final projected gradient 0.0598606
## final function value -1635.49
## 
## F = -1635.49
## Warning:  more than 10 function and gradient evaluations
##    in the last line search
## final  value -1635.491845 
## converged
## Optimisation using starting value 11/12
## N = 6, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate     0  f=      -1370.2  |proj g|=           20
## At iterate    10  f =      -1635.4  |proj g|=        1.0976
## 
## iterations 15
## function evaluations 21
## segments explored during Cauchy searches 20
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 0
## norm of the final projected gradient 0.00186049
## final function value -1635.49
## 
## F = -1635.49
## final  value -1635.491833 
## converged
## Optimisation using starting value 12/12
## N = 6, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate     0  f=      -1373.5  |proj g|=           20
## At iterate    10  f =      -1635.5  |proj g|=       0.19953
## 
## iterations 15
## function evaluations 31
## segments explored during Cauchy searches 19
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 0
## norm of the final projected gradient 0.0461308
## final function value -1635.49
## 
## F = -1635.49
## Warning:  more than 10 function and gradient evaluations
##    in the last line search
## final  value -1635.491840 
## converged
print(est.denver.model.B)
## Optimisation for STmodel with 12 starting points.
##   Results: 4 converged, 8 not converged, 0 failed.
##   Best result for starting point 4, optimisation has converged
## 
## No fixed parameters.
## 
## Estimated parameters for all starting point(s):
##                                         [,1]            [,2]             [,3]
## gamma.bc_st_no2                 0.0060118836   0.00596005607   0.005955238637
## gamma.bc_st_smk                 0.0521824404   0.05195395316   0.051963119374
## alpha.const.(Intercept)         0.1912785980   0.25013512306   0.250232500893
## alpha.const.tree_cover_100      0.0032725799   0.00858563403   0.008600000285
## alpha.const.impervious_2500     0.0086346267   0.01913402640   0.019148742535
## alpha.const.low_int_100         0.0093760168   0.01430191351   0.014301194760
## alpha.const.med_int_50         -0.0039613701  -0.00421505227  -0.004210495232
## alpha.const.med_int_2500        0.0171963057  -0.00094776366  -0.000943197470
## alpha.const.high_int_100        0.0014554826   0.01229113752   0.012293428681
## alpha.const.pop_den_50         -0.0158271837  -0.01315055513  -0.013153454350
## alpha.const.dist_m_cafos       -0.0196381194  -0.02064697856  -0.020650426238
## alpha.const.dist_m_og_wells    -0.0014823607   0.00410607533   0.004094448529
## alpha.const.dist_m_wwtp         0.0061514019   0.00280361913   0.002820074602
## alpha.const.aadt_100            0.0344456466   0.01913393117   0.019130992908
## alpha.const.aadt_250            0.0216541261   0.02327931343   0.023287934150
## alpha.const.aadt_1000          -0.0082169471  -0.00884375176  -0.008843479197
## alpha.V1.(Intercept)            0.1640948824   0.18448242880   0.184486002539
## alpha.V1.tree_cover_100        -0.0149371422  -0.01216815574  -0.012157621514
## alpha.V1.impervious_2500        0.0046700767  -0.00589806167  -0.005876364566
## alpha.V1.low_int_100           -0.0099590774  -0.00751872533  -0.007518516324
## alpha.V1.med_int_50            -0.0145381299  -0.01258187668  -0.012580731420
## alpha.V1.med_int_2500           0.0007319780   0.00362165651   0.003607057482
## alpha.V1.high_int_100           0.0193619626   0.01805976475   0.018074702267
## alpha.V1.pop_den_50             0.0027890733   0.00271971642   0.002727439623
## alpha.V1.dist_m_cafos           0.0059170892  -0.00001866864  -0.000001210275
## alpha.V1.dist_m_og_wells       -0.0016883212   0.00379434591   0.003797105516
## alpha.V1.dist_m_wwtp           -0.0027537848  -0.00835900122  -0.008348554976
## alpha.V1.aadt_100              -0.0244963327  -0.03129890664  -0.031308202681
## alpha.V1.aadt_250              -0.0206482343  -0.00905272933  -0.009062198058
## alpha.V1.aadt_1000              0.0044713009   0.00044568561   0.000442089992
## alpha.V2.(Intercept)            0.1604857126   0.18506394002   0.185100939006
## alpha.V2.tree_cover_100        -0.0038598984  -0.00140760762  -0.001396523996
## alpha.V2.impervious_2500       -0.0292271566  -0.02688644737  -0.026870279271
## alpha.V2.low_int_100           -0.0008278762   0.00081574734   0.000813085668
## alpha.V2.med_int_50            -0.0085606701  -0.00852572347  -0.008520145449
## alpha.V2.med_int_2500           0.0213995519   0.01609592640   0.016090534357
## alpha.V2.high_int_100          -0.0115300275  -0.00531046255  -0.005307474586
## alpha.V2.pop_den_50             0.0070030441   0.00627161540   0.006274912103
## alpha.V2.dist_m_cafos           0.0143984442   0.01517681980   0.015169579556
## alpha.V2.dist_m_og_wells       -0.0018465363   0.00316452234   0.003156172950
## alpha.V2.dist_m_wwtp            0.0009020415  -0.00309514261  -0.003081324416
## alpha.V2.aadt_100              -0.0189470731  -0.02803138572  -0.028037993207
## alpha.V2.aadt_250               0.0107539488   0.01242409203   0.012422670722
## alpha.V2.aadt_1000             -0.0056846619  -0.00812879018  -0.008124952189
## log.nugget.const.iid          -14.2075145954 -14.99853953428 -14.998066086823
## log.nugget.V1.iid             -13.3885208673 -13.25977435576 -13.811919503385
## log.nugget.V2.iid             -15.0000000000 -14.99612179953 -14.990498086304
## nu.log.range.exp                0.0000000000  13.75456682509  13.762410630096
## nu.log.sill.exp                -4.6993271432  -3.00667112679  -3.006919056540
## nu.log.nugget.(Intercept).exp  -4.2038352592  -4.76493602901  -4.764763069367
##                                         [,4]           [,5]          [,6]
## gamma.bc_st_no2                 0.0062756184   0.0062750888  0.0062448184
## gamma.bc_st_smk                 0.0518645345   0.0518633379  0.0518689417
## alpha.const.(Intercept)         0.2490535778   0.2490551735  0.2494632803
## alpha.const.tree_cover_100      0.0048752209   0.0048744305  0.0049646999
## alpha.const.impervious_2500     0.0196052934   0.0196284557  0.0197915602
## alpha.const.low_int_100         0.0139807586   0.0139840762  0.0140120667
## alpha.const.med_int_50         -0.0090675828  -0.0090833998 -0.0091188953
## alpha.const.med_int_2500       -0.0042252843  -0.0042602214 -0.0044315319
## alpha.const.high_int_100        0.0098436326   0.0098418272  0.0098587736
## alpha.const.pop_den_50         -0.0097277268  -0.0097230856 -0.0097633578
## alpha.const.dist_m_cafos       -0.0212607744  -0.0212768491 -0.0213592746
## alpha.const.dist_m_og_wells     0.0063053930   0.0063061812  0.0062501046
## alpha.const.dist_m_wwtp        -0.0014448112  -0.0014503791 -0.0014192174
## alpha.const.aadt_100            0.0198837181   0.0198758207  0.0198041276
## alpha.const.aadt_250            0.0232382337   0.0232435263  0.0233246181
## alpha.const.aadt_1000          -0.0102055575  -0.0102042173 -0.0101736236
## alpha.V1.(Intercept)            0.1796396617   0.1796034928  0.1797166847
## alpha.V1.tree_cover_100        -0.0095807532  -0.0095857651 -0.0095930706
## alpha.V1.impervious_2500       -0.0026502125  -0.0026663502 -0.0026861356
## alpha.V1.low_int_100           -0.0063794936  -0.0063868015 -0.0064147407
## alpha.V1.med_int_50            -0.0115235700  -0.0115134344 -0.0114282183
## alpha.V1.med_int_2500          -0.0001957307  -0.0001779400 -0.0000751741
## alpha.V1.high_int_100           0.0210467647   0.0210293676  0.0209677001
## alpha.V1.pop_den_50             0.0030558366   0.0030481096  0.0029986132
## alpha.V1.dist_m_cafos          -0.0014346281  -0.0014369914 -0.0014559278
## alpha.V1.dist_m_og_wells        0.0011285833   0.0011302192  0.0012182994
## alpha.V1.dist_m_wwtp           -0.0045033269  -0.0045077458 -0.0045269901
## alpha.V1.aadt_100              -0.0254282735  -0.0254231814 -0.0254880303
## alpha.V1.aadt_250              -0.0107957045  -0.0107872255 -0.0107660633
## alpha.V1.aadt_1000             -0.0003513547  -0.0003443569 -0.0003134788
## alpha.V2.(Intercept)            0.1837226965   0.1837249799  0.1838253301
## alpha.V2.tree_cover_100        -0.0032709767  -0.0032725068 -0.0032224787
## alpha.V2.impervious_2500       -0.0224524331  -0.0224434950 -0.0223627890
## alpha.V2.low_int_100            0.0020186127   0.0020206665  0.0020292072
## alpha.V2.med_int_50            -0.0113583401  -0.0113640954 -0.0113460420
## alpha.V2.med_int_2500           0.0096917376   0.0096847193  0.0096683086
## alpha.V2.high_int_100          -0.0053929727  -0.0053929997 -0.0053821449
## alpha.V2.pop_den_50             0.0087220406   0.0087151215  0.0086414050
## alpha.V2.dist_m_cafos           0.0125478935   0.0125360189  0.0124589374
## alpha.V2.dist_m_og_wells        0.0035149242   0.0035251114  0.0035775837
## alpha.V2.dist_m_wwtp           -0.0032472277  -0.0032520674 -0.0032526984
## alpha.V2.aadt_100              -0.0273235923  -0.0273252702 -0.0273759384
## alpha.V2.aadt_250               0.0125332866   0.0125343813  0.0125506316
## alpha.V2.aadt_1000             -0.0082224575  -0.0082276268 -0.0082511959
## log.nugget.const.iid          -10.4632900888 -10.3985631570 -9.9979407159
## log.nugget.V1.iid              -7.0285457395  -7.0243931348 -7.0258773232
## log.nugget.V2.iid              -7.7309663192  -7.7317201739 -7.7474076123
## nu.log.range.exp               13.5125634739  13.5098031295 13.5109546812
## nu.log.sill.exp                -2.9784990884  -2.9783780489 -2.9780490360
## nu.log.nugget.(Intercept).exp  -4.9340326598  -4.9344286674 -4.9352276245
##                                         [,7]            [,8]           [,9]
## gamma.bc_st_no2                 0.0062750305   0.00624754768  0.00624517200
## gamma.bc_st_smk                 0.0518629046   0.05186815364  0.05186799148
## alpha.const.(Intercept)         0.2490496581   0.24942294174  0.24945251190
## alpha.const.tree_cover_100      0.0048788656   0.00495761154  0.00496455640
## alpha.const.impervious_2500     0.0196042830   0.01977807254  0.01979078798
## alpha.const.low_int_100         0.0139810141   0.01400994065  0.01401226517
## alpha.const.med_int_50         -0.0090654055  -0.00911525896 -0.00911852466
## alpha.const.med_int_2500       -0.0042249847  -0.00441729279 -0.00443215946
## alpha.const.high_int_100        0.0098442364   0.00985798842  0.00985905250
## alpha.const.pop_den_50         -0.0097324016  -0.00976044755 -0.00976428759
## alpha.const.dist_m_cafos       -0.0212606928  -0.02135274314 -0.02135967695
## alpha.const.dist_m_og_wells     0.0063029166   0.00625448630  0.00625008615
## alpha.const.dist_m_wwtp        -0.0014443857  -0.00142087749 -0.00141947406
## alpha.const.aadt_100            0.0198820425   0.01981001382  0.01980396759
## alpha.const.aadt_250            0.0232404112   0.02331764231  0.02332424882
## alpha.const.aadt_1000          -0.0102036016  -0.01017580725 -0.01017292717
## alpha.V1.(Intercept)            0.1796525981   0.17970823595  0.17971635560
## alpha.V1.tree_cover_100        -0.0095817610  -0.00959386626 -0.00959521360
## alpha.V1.impervious_2500       -0.0026517806  -0.00268701723 -0.00269063387
## alpha.V1.low_int_100           -0.0063784085  -0.00641321983 -0.00641533893
## alpha.V1.med_int_50            -0.0115203868  -0.01143558341 -0.01142777810
## alpha.V1.med_int_2500          -0.0001905905  -0.00008132414 -0.00007091681
## alpha.V1.high_int_100           0.0210459615   0.02097097055  0.02096485396
## alpha.V1.pop_den_50             0.0030534950   0.00300235265  0.00299744584
## alpha.V1.dist_m_cafos          -0.0014387133  -0.00145376537 -0.00145682140
## alpha.V1.dist_m_og_wells        0.0011348363   0.00121167969  0.00121976868
## alpha.V1.dist_m_wwtp           -0.0045044776  -0.00452732568 -0.00452936877
## alpha.V1.aadt_100              -0.0254339827  -0.02548402101 -0.02548962933
## alpha.V1.aadt_250              -0.0107952680  -0.01076688982 -0.01076441074
## alpha.V1.aadt_1000             -0.0003501826  -0.00031577454 -0.00031259535
## alpha.V2.(Intercept)            0.1837201618   0.18381566669  0.18382257167
## alpha.V2.tree_cover_100        -0.0032693538  -0.00322654205 -0.00322298325
## alpha.V2.impervious_2500       -0.0224532625  -0.02237127152 -0.02236530503
## alpha.V2.low_int_100            0.0020185719   0.00202819139  0.00202903497
## alpha.V2.med_int_50            -0.0113557183  -0.01134716779 -0.01134570891
## alpha.V2.med_int_2500           0.0096941148   0.00967168608  0.00967055008
## alpha.V2.high_int_100          -0.0053933123  -0.00538316208 -0.00538260579
## alpha.V2.pop_den_50             0.0087181759   0.00864682945  0.00863992996
## alpha.V2.dist_m_cafos           0.0125471507   0.01246616230  0.01245955972
## alpha.V2.dist_m_og_wells        0.0035164838   0.00357346862  0.00357858028
## alpha.V2.dist_m_wwtp           -0.0032481932  -0.00325304734 -0.00325403892
## alpha.V2.aadt_100              -0.0273253755  -0.02737162877 -0.02737577550
## alpha.V2.aadt_250               0.0125341258   0.01254927405  0.01255077775
## alpha.V2.aadt_1000             -0.0082225089  -0.00824950170 -0.00825156424
## log.nugget.const.iid          -10.4481949545 -10.02637857089 -9.99617297360
## log.nugget.V1.iid              -7.0300419363  -7.02605082295 -7.02622038251
## log.nugget.V2.iid              -7.7315355639  -7.74665482316 -7.74794236783
## nu.log.range.exp               13.5121556541  13.51065516608 13.51030811048
## nu.log.sill.exp                -2.9784997305  -2.97815642597 -2.97807591297
## nu.log.nugget.(Intercept).exp  -4.9340481552  -4.93516126191 -4.93527648027
##                                        [,10]          [,11]          [,12]
## gamma.bc_st_no2                 0.0063171711   0.0063168583   0.0063113077
## gamma.bc_st_smk                 0.0518590513   0.0518631206   0.0518749984
## alpha.const.(Intercept)         0.2485173757   0.2485413475   0.2486812126
## alpha.const.tree_cover_100      0.0047559128   0.0047561336   0.0047672588
## alpha.const.impervious_2500     0.0193399493   0.0193531093   0.0193839872
## alpha.const.low_int_100         0.0139381060   0.0139393168   0.0139435416
## alpha.const.med_int_50         -0.0089886588  -0.0089936373  -0.0089963190
## alpha.const.med_int_2500       -0.0039321095  -0.0039414894  -0.0039561026
## alpha.const.high_int_100        0.0098228588   0.0098238091   0.0098299831
## alpha.const.pop_den_50         -0.0096881146  -0.0096810069  -0.0096742703
## alpha.const.dist_m_cafos       -0.0211182900  -0.0211236440  -0.0211321894
## alpha.const.dist_m_og_wells     0.0063769747   0.0063774050   0.0063714233
## alpha.const.dist_m_wwtp        -0.0014812408  -0.0014771607  -0.0014596871
## alpha.const.aadt_100            0.0199933751   0.0199919591   0.0199874581
## alpha.const.aadt_250            0.0231216219   0.0231228210   0.0231301077
## alpha.const.aadt_1000          -0.0102441085  -0.0102462022  -0.0102482989
## alpha.V1.(Intercept)            0.1795708992   0.1795570796   0.1795671190
## alpha.V1.tree_cover_100        -0.0095625102  -0.0095602429  -0.0095541761
## alpha.V1.impervious_2500       -0.0025956207  -0.0025908271  -0.0025747943
## alpha.V1.low_int_100           -0.0063244373  -0.0063289545  -0.0063375917
## alpha.V1.med_int_50            -0.0116509561  -0.0116532175  -0.0116568871
## alpha.V1.med_int_2500          -0.0003597753  -0.0003655688  -0.0003772143
## alpha.V1.high_int_100           0.0211614029   0.0211601862   0.0211612946
## alpha.V1.pop_den_50             0.0031321987   0.0031348258   0.0031406720
## alpha.V1.dist_m_cafos          -0.0014137351  -0.0014040404  -0.0013821293
## alpha.V1.dist_m_og_wells        0.0010156107   0.0010083986   0.0010010728
## alpha.V1.dist_m_wwtp           -0.0044677063  -0.0044667258  -0.0044660557
## alpha.V1.aadt_100              -0.0253585159  -0.0253526404  -0.0253508028
## alpha.V1.aadt_250              -0.0108364819  -0.0108368589  -0.0108408111
## alpha.V1.aadt_1000             -0.0004058246  -0.0004074424  -0.0004118690
## alpha.V2.(Intercept)            0.1835930005   0.1836037312   0.1836445815
## alpha.V2.tree_cover_100        -0.0033363676  -0.0033349391  -0.0033253831
## alpha.V2.impervious_2500       -0.0225745467  -0.0225672019  -0.0225518680
## alpha.V2.low_int_100            0.0020053069   0.0020054011   0.0020050406
## alpha.V2.med_int_50            -0.0113689933  -0.0113722966  -0.0113732523
## alpha.V2.med_int_2500           0.0097242073   0.0097193770   0.0097146054
## alpha.V2.high_int_100          -0.0054088191  -0.0054068723  -0.0054014877
## alpha.V2.pop_den_50             0.0088287099   0.0088324299   0.0088371225
## alpha.V2.dist_m_cafos           0.0126710514   0.0126685543   0.0126643029
## alpha.V2.dist_m_og_wells        0.0034284932   0.0034273196   0.0034236107
## alpha.V2.dist_m_wwtp           -0.0032402629  -0.0032366303  -0.0032262324
## alpha.V2.aadt_100              -0.0272558331  -0.0272558188  -0.0272598777
## alpha.V2.aadt_250               0.0125113744   0.0125105768   0.0125097559
## alpha.V2.aadt_1000             -0.0081800766  -0.0081813330  -0.0081832037
## log.nugget.const.iid          -11.9828653970 -11.9963559520 -12.0001162114
## log.nugget.V1.iid              -7.0355605905  -7.0331561647  -7.0306406707
## log.nugget.V2.iid              -7.7081729398  -7.7085438681  -7.7106093212
## nu.log.range.exp               13.5151755391  13.5167035286  13.5224235635
## nu.log.sill.exp                -2.9785879051  -2.9790274847  -2.9789261257
## nu.log.nugget.(Intercept).exp  -4.9320850646  -4.9321253281  -4.9314759991
## 
## Function value(s):
##  [1] 1394.207 1622.137 1622.131 1635.494 1635.494 1635.492 1635.494 1635.493
##  [9] 1635.492 1635.492 1635.492 1635.492

3.2.3 Cross-validation

Define the CV groups

set.seed(1000)

site_idsB <- colnames(denver.model.B$D.beta)[which(str_detect(colnames(denver.model.B$D.beta), "d_") == T)]
site_idsB
##  [1] "d_1"  "d_10" "d_11" "d_12" "d_13" "d_14" "d_15" "d_16" "d_17" "d_18"
## [11] "d_19" "d_2"  "d_20" "d_21" "d_22" "d_23" "d_25" "d_26" "d_28" "d_29"
## [21] "d_3"  "d_31" "d_32" "d_33" "d_34" "d_35" "d_36" "d_37" "d_38" "d_39"
## [31] "d_4"  "d_40" "d_41" "d_42" "d_43" "d_44" "d_45" "d_46" "d_47" "d_48"
## [41] "d_49" "d_5"  "d_50" "d_51" "d_52" "d_53" "d_54" "d_55" "d_56" "d_57"
## [51] "d_58" "d_6"  "d_60" "d_61" "d_62" "d_63" "d_64" "d_65" "d_66" "d_67"
## [61] "d_68" "d_7"  "d_8"
Ind.cv.B <- createCV(denver.model.B, groups = 10, #min.dist = .1,
                     subset = site_idsB)

ID.cv.B <- sapply(split(denver.model.B$obs$ID, Ind.cv.B), unique)
print(sapply(ID.cv.B, length))
##  0  1  2  3  4  5  6  7  8  9 10 
##  1  7  7  7  6  6  6  6  6  6  6
table(Ind.cv.B)
## Ind.cv.B
##   0   1   2   3   4   5   6   7   8   9  10 
## 231  75  80  97  73  64  93  80  81  66  81
I.col.B <- apply(sapply(ID.cv.B, function(x) denver.model.B$locations$ID%in% x), 1,
                 function(x) if(sum(x)==1) which(x) else 0)
names(I.col.B) <- denver.model.B$locations$ID
print(I.col.B)
## central     d_1    d_10    d_11    d_12    d_13    d_14    d_15    d_16    d_17 
##       1       6       2       3       4       4       4       9       4      10 
##    d_18    d_19     d_2    d_20    d_21    d_22    d_23    d_25    d_26    d_28 
##       8       5      11       3       2      11       4       9      11       9 
##    d_29     d_3    d_31    d_32    d_33    d_34    d_35    d_36    d_37    d_38 
##       2       4       6       8       5       5      10       8       3       8 
##    d_39     d_4    d_40    d_41    d_42    d_43    d_44    d_45    d_46    d_47 
##       6       7       7       3       7       8      11       3       9       5 
##    d_48    d_49     d_5    d_50    d_51    d_52    d_53    d_54    d_55    d_56 
##       9       4       6       9      10       3       7       5       7       6 
##    d_57    d_58     d_6    d_60    d_61    d_62    d_63    d_64    d_65    d_66 
##       2       2      11       5       2      11      10       3       8       6 
##    d_67    d_68     d_7     d_8    d_24    d_27    d_30    d_59     d_9 
##       7      10       2      10       0       0       0       0       0
par(mfrow=c(1,1))
plot(denver.model.B$locations$long,
     denver.model.B$locations$lat,
     pch=23+floor(I.col.B/max(I.col.B)+.5), bg=I.col.B,
     xlab="Longitude", ylab="Latitude")

map("county", "colorado", col="#FFFF0055",fill=TRUE, add=TRUE)
## [[1]]
## NULL

ID starting conditions using previous model without CV:

x.init.B.cv <- coef(est.denver.model.B, pars="cov")[,c("par","init")]
x.init.B.cv

Run the model with cross validation

est.denver.B.cv <- estimateCV(denver.model.B, x.init.B.cv, Ind.cv.B)
## 
## ***************************
## Estimation of CV-set 1/10
## Optimisation using starting value 1/2
## N = 6, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate     0  f=      -1484.4  |proj g|=       13.035
## At iterate    10  f =      -1485.3  |proj g|=      0.076883
## At iterate    20  f =      -1485.3  |proj g|=     0.0093172
## 
## iterations 21
## function evaluations 35
## segments explored during Cauchy searches 21
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 0
## norm of the final projected gradient 0.00931736
## final function value -1485.34
## 
## F = -1485.34
## final  value -1485.339955 
## converged
## 
## Optimisation using starting value 2/2
## N = 6, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate     0  f=      -1105.4  |proj g|=           12
## At iterate    10  f =      -1479.4  |proj g|=       0.73963
## At iterate    20  f =      -1483.4  |proj g|=        13.853
## At iterate    30  f =      -1485.3  |proj g|=       0.94886
## 
## iterations 39
## function evaluations 64
## segments explored during Cauchy searches 42
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 0
## norm of the final projected gradient 0.00705453
## final function value -1485.34
## 
## F = -1485.34
## Warning:  more than 10 function and gradient evaluations
##    in the last line search
## final  value -1485.343747 
## converged
## 
## 
## ***************************
## Estimation of CV-set 2/10
## Optimisation using starting value 1/2
## N = 6, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate     0  f=        -1500  |proj g|=       8.8707
## At iterate    10  f =      -1500.3  |proj g|=       0.39016
## At iterate    20  f =      -1500.3  |proj g|=       0.24659
## At iterate    30  f =      -1500.4  |proj g|=        0.5606
## At iterate    40  f =      -1500.4  |proj g|=      0.071363
## 
## iterations 45
## function evaluations 54
## segments explored during Cauchy searches 46
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 1
## norm of the final projected gradient 0.00240616
## final function value -1500.38
## 
## F = -1500.38
## final  value -1500.383951 
## converged
## 
## Optimisation using starting value 2/2
## N = 6, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate     0  f=      -1104.7  |proj g|=           12
## At iterate    10  f =      -1424.8  |proj g|=        11.786
## At iterate    20  f =      -1487.1  |proj g|=       0.93927
## At iterate    30  f =      -1489.7  |proj g|=        2.6261
## At iterate    40  f =      -1490.3  |proj g|=        1.4331
## At iterate    50  f =        -1494  |proj g|=        1.6876
## At iterate    60  f =      -1494.5  |proj g|=       0.13324
## At iterate    70  f =      -1494.5  |proj g|=       0.10735
## At iterate    80  f =      -1499.1  |proj g|=        4.7371
## At iterate    90  f =      -1500.3  |proj g|=       0.79519
## At iterate   100  f =      -1500.4  |proj g|=        0.0523
## 
## iterations 104
## function evaluations 137
## segments explored during Cauchy searches 107
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 0
## norm of the final projected gradient 0.143867
## final function value -1500.38
## 
## F = -1500.38
## final  value -1500.382816 
## converged
## 
## 
## ***************************
## Estimation of CV-set 3/10
## Optimisation using starting value 1/2
## N = 6, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate     0  f=      -1467.3  |proj g|=       10.066
## At iterate    10  f =      -1467.6  |proj g|=       0.15397
## 
## iterations 14
## function evaluations 18
## segments explored during Cauchy searches 15
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 0
## norm of the final projected gradient 0.0369096
## final function value -1467.62
## 
## F = -1467.62
## final  value -1467.616601 
## converged
## 
## Optimisation using starting value 2/2
## N = 6, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate     0  f=      -1089.1  |proj g|=           12
## At iterate    10  f =        -1400  |proj g|=        12.375
## At iterate    20  f =        -1458  |proj g|=        2.2193
## At iterate    30  f =      -1460.4  |proj g|=       0.14393
## At iterate    40  f =      -1463.9  |proj g|=        1.0706
## At iterate    50  f =      -1463.9  |proj g|=      0.013697
## At iterate    60  f =      -1466.9  |proj g|=        1.7069
## At iterate    70  f =      -1467.6  |proj g|=       0.45824
## At iterate    80  f =      -1467.6  |proj g|=        0.1933
## 
## iterations 82
## function evaluations 122
## segments explored during Cauchy searches 85
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 0
## norm of the final projected gradient 0.0399222
## final function value -1467.62
## 
## F = -1467.62
## final  value -1467.618744 
## converged
## 
## 
## ***************************
## Estimation of CV-set 4/10
## Optimisation using starting value 1/2
## N = 6, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate     0  f=      -1481.7  |proj g|=       17.346
## At iterate    10  f =      -1482.6  |proj g|=        0.3585
## At iterate    20  f =      -1482.6  |proj g|=      0.032367
## At iterate    30  f =      -1482.7  |proj g|=       0.36671
## At iterate    40  f =      -1482.7  |proj g|=       0.12535
## At iterate    50  f =      -1482.7  |proj g|=       0.39405
## 
## iterations 56
## function evaluations 78
## segments explored during Cauchy searches 56
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 0
## norm of the final projected gradient 0.0209832
## final function value -1482.67
## 
## F = -1482.67
## Warning:  more than 10 function and gradient evaluations
##    in the last line search
## final  value -1482.672168 
## converged
## 
## Optimisation using starting value 2/2
## N = 6, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate     0  f=      -1107.5  |proj g|=           12
## At iterate    10  f =      -1416.2  |proj g|=        11.607
## At iterate    20  f =      -1473.2  |proj g|=        2.2818
## At iterate    30  f =      -1475.1  |proj g|=        9.4907
## At iterate    40  f =      -1475.7  |proj g|=       0.02604
## At iterate    50  f =      -1477.7  |proj g|=        1.0248
## At iterate    60  f =      -1477.8  |proj g|=       0.07492
## At iterate    70  f =      -1482.2  |proj g|=        2.4183
## At iterate    80  f =      -1482.5  |proj g|=       0.80487
## At iterate    90  f =      -1482.6  |proj g|=       0.37974
## At iterate   100  f =      -1482.7  |proj g|=       0.38035
## At iterate   110  f =      -1482.7  |proj g|=      0.057702
## 
## iterations 119
## function evaluations 168
## segments explored during Cauchy searches 122
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 1
## norm of the final projected gradient 0.0145089
## final function value -1482.67
## 
## F = -1482.67
## final  value -1482.672380 
## converged
## 
## 
## ***************************
## Estimation of CV-set 5/10
## Optimisation using starting value 1/2
## N = 6, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate     0  f=        -1506  |proj g|=       10.656
## At iterate    10  f =      -1506.5  |proj g|=      0.059847
## At iterate    20  f =      -1506.6  |proj g|=        0.3413
## At iterate    30  f =      -1506.6  |proj g|=       0.15811
## Bad direction in the line search;
##    refresh the lbfgs memory and restart the iteration.
## Line search cannot locate an adequate point after 20 function
## and gradient evaluations
## final  value -1506.581892 
## stopped after 35 iterations
## 
## Optimisation using starting value 2/2
## N = 6, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate     0  f=      -1120.3  |proj g|=           12
## At iterate    10  f =      -1493.7  |proj g|=        12.155
## At iterate    20  f =      -1499.8  |proj g|=         4.429
## At iterate    30  f =      -1506.3  |proj g|=        1.8316
## At iterate    40  f =      -1506.6  |proj g|=       0.18166
## At iterate    50  f =      -1506.6  |proj g|=       0.39657
## At iterate    60  f =      -1506.6  |proj g|=      0.098565
## 
## iterations 65
## function evaluations 97
## segments explored during Cauchy searches 68
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 1
## norm of the final projected gradient 0.102392
## final function value -1506.58
## 
## F = -1506.58
## Warning:  more than 10 function and gradient evaluations
##    in the last line search
## final  value -1506.582174 
## converged
## 
## 
## ***************************
## Estimation of CV-set 6/10
## Optimisation using starting value 1/2
## N = 6, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate     0  f=      -1481.1  |proj g|=       10.066
## At iterate    10  f =      -1481.5  |proj g|=       0.06454
## At iterate    20  f =      -1481.6  |proj g|=       0.90862
## At iterate    30  f =      -1481.6  |proj g|=       0.07344
## Bad direction in the line search;
##    refresh the lbfgs memory and restart the iteration.
## Line search cannot locate an adequate point after 20 function
## and gradient evaluations
## final  value -1481.582531 
## stopped after 35 iterations
## 
## Optimisation using starting value 2/2
## N = 6, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate     0  f=      -1086.4  |proj g|=           12
## At iterate    10  f =      -1471.7  |proj g|=        3.8336
## At iterate    20  f =      -1479.3  |proj g|=        2.7325
## At iterate    30  f =      -1481.3  |proj g|=        5.9125
## At iterate    40  f =      -1481.5  |proj g|=       0.33309
## At iterate    50  f =      -1481.6  |proj g|=       0.25228
## At iterate    60  f =      -1481.6  |proj g|=      0.034288
## At iterate    70  f =      -1481.6  |proj g|=      0.027797
## 
## iterations 74
## function evaluations 90
## segments explored during Cauchy searches 77
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 1
## norm of the final projected gradient 0.0336023
## final function value -1481.58
## 
## F = -1481.58
## final  value -1481.582484 
## converged
## 
## 
## ***************************
## Estimation of CV-set 7/10
## Optimisation using starting value 1/2
## N = 6, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate     0  f=      -1491.1  |proj g|=       2.1306
## At iterate    10  f =      -1492.4  |proj g|=       0.95479
## At iterate    20  f =      -1492.4  |proj g|=       0.12235
## At iterate    30  f =      -1492.5  |proj g|=       0.46536
## At iterate    40  f =      -1492.5  |proj g|=       0.33101
## At iterate    50  f =      -1492.5  |proj g|=     0.0067843
## 
## iterations 50
## function evaluations 57
## segments explored during Cauchy searches 51
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 1
## norm of the final projected gradient 0.00678434
## final function value -1492.55
## 
## F = -1492.55
## final  value -1492.549537 
## converged
## 
## Optimisation using starting value 2/2
## N = 6, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate     0  f=      -1099.7  |proj g|=           12
## At iterate    10  f =      -1484.8  |proj g|=          1.04
## At iterate    20  f =      -1486.8  |proj g|=        7.7877
## At iterate    30  f =      -1489.2  |proj g|=       0.57899
## At iterate    40  f =      -1492.4  |proj g|=       0.33942
## At iterate    50  f =      -1492.5  |proj g|=        0.3765
## 
## iterations 59
## function evaluations 90
## segments explored during Cauchy searches 62
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 1
## norm of the final projected gradient 0.0306232
## final function value -1492.55
## 
## F = -1492.55
## final  value -1492.549474 
## converged
## 
## 
## ***************************
## Estimation of CV-set 8/10
## Optimisation using starting value 1/2
## N = 6, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate     0  f=      -1479.5  |proj g|=       9.3568
## At iterate    10  f =      -1480.6  |proj g|=        1.2495
## At iterate    20  f =      -1480.7  |proj g|=       0.48288
## At iterate    30  f =        -1481  |proj g|=         0.188
## At iterate    40  f =        -1481  |proj g|=     0.0058638
## 
## iterations 41
## function evaluations 61
## segments explored during Cauchy searches 41
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 0
## norm of the final projected gradient 0.00586381
## final function value -1481
## 
## F = -1481
## Warning:  more than 10 function and gradient evaluations
##    in the last line search
## final  value -1481.002261 
## converged
## 
## Optimisation using starting value 2/2
## N = 6, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate     0  f=      -1097.9  |proj g|=           12
## At iterate    10  f =      -1473.3  |proj g|=        12.094
## At iterate    20  f =      -1478.8  |proj g|=        2.3195
## At iterate    30  f =        -1480  |proj g|=       0.67451
## At iterate    40  f =        -1481  |proj g|=       0.15932
## 
## iterations 49
## function evaluations 77
## segments explored during Cauchy searches 52
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 0
## norm of the final projected gradient 0.016169
## final function value -1481
## 
## F = -1481
## Warning:  more than 10 function and gradient evaluations
##    in the last line search
## final  value -1481.002269 
## converged
## 
## 
## ***************************
## Estimation of CV-set 9/10
## Optimisation using starting value 1/2
## N = 6, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate     0  f=      -1500.8  |proj g|=       11.701
## At iterate    10  f =      -1501.3  |proj g|=       0.13133
## 
## iterations 15
## function evaluations 31
## segments explored during Cauchy searches 15
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 0
## norm of the final projected gradient 0.10318
## final function value -1501.33
## 
## F = -1501.33
## Warning:  more than 10 function and gradient evaluations
##    in the last line search
## final  value -1501.329327 
## converged
## 
## Optimisation using starting value 2/2
## N = 6, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate     0  f=      -1117.3  |proj g|=           12
## At iterate    10  f =      -1493.4  |proj g|=        6.7681
## At iterate    20  f =      -1497.9  |proj g|=         14.57
## At iterate    30  f =      -1501.2  |proj g|=       0.35627
## At iterate    40  f =      -1501.3  |proj g|=       0.62611
## At iterate    50  f =      -1501.3  |proj g|=      0.036213
## 
## iterations 56
## function evaluations 80
## segments explored during Cauchy searches 59
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 0
## norm of the final projected gradient 0.0143426
## final function value -1501.33
## 
## F = -1501.33
## final  value -1501.334466 
## converged
## 
## 
## ***************************
## Estimation of CV-set 10/10
## Optimisation using starting value 1/2
## N = 6, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate     0  f=      -1499.1  |proj g|=       10.066
## At iterate    10  f =      -1500.3  |proj g|=        0.2129
## At iterate    20  f =      -1500.3  |proj g|=      0.059763
## At iterate    30  f =      -1500.3  |proj g|=       0.20161
## At iterate    40  f =      -1500.3  |proj g|=      0.058964
## At iterate    50  f =      -1500.3  |proj g|=       0.00997
## 
## iterations 52
## function evaluations 68
## segments explored during Cauchy searches 53
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 1
## norm of the final projected gradient 0.0152949
## final function value -1500.3
## 
## F = -1500.3
## Warning:  more than 10 function and gradient evaluations
##    in the last line search
## final  value -1500.301602 
## converged
## 
## Optimisation using starting value 2/2
## N = 6, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate     0  f=      -1099.5  |proj g|=           12
## At iterate    10  f =      -1488.4  |proj g|=        5.2544
## At iterate    20  f =      -1491.7  |proj g|=        5.0281
## At iterate    30  f =      -1493.6  |proj g|=        10.355
## At iterate    40  f =        -1498  |proj g|=        1.1003
## At iterate    50  f =        -1500  |proj g|=       0.78914
## At iterate    60  f =      -1500.3  |proj g|=        0.9166
## At iterate    70  f =      -1500.3  |proj g|=       0.31545
## At iterate    80  f =      -1500.3  |proj g|=      0.043618
## 
## iterations 83
## function evaluations 128
## segments explored during Cauchy searches 86
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 0
## norm of the final projected gradient 0.0248369
## final function value -1500.3
## 
## F = -1500.3
## final  value -1500.301353 
## converged
## 
print(est.denver.B.cv)
## Cross-validation parameter estimation for STmodel
##   with 10 CV-groups and 2 starting points.
##   Results: 7 converged, 3 not converged.
## 
## No fixed parameters.
## 
## Estimated function values and convergence info:
##       value convergence  conv     eigen.min eigen.all.min
## 1  1485.344        TRUE  TRUE  0.0036525174            NA
## 2  1500.384        TRUE FALSE -0.0022060915            NA
## 3  1467.619        TRUE  TRUE  0.0028793325            NA
## 4  1482.672        TRUE FALSE -0.0001855608            NA
## 5  1506.582        TRUE  TRUE  0.0002412768            NA
## 6  1481.582        TRUE  TRUE  0.0001981437            NA
## 7  1492.550        TRUE  TRUE  0.0009942168            NA
## 8  1481.002        TRUE  TRUE  0.5570127079            NA
## 9  1501.334        TRUE  TRUE  0.0439724903            NA
## 10 1500.302        TRUE FALSE -0.0008568626            NA
par(mfrow=c(1,1), mar=c(13.5,2.5,.5,.5), las=2)
with(coef(est.denver.model.B, pars="all"),
     plotCI((1:length(par))+.3, par, uiw=1.96*sd,
            col=2, xlab="", xaxt="n", ylab=""))
boxplot(est.denver.B.cv, "all", boxwex=.4, col="grey", add=TRUE)

#' Save the results as an .rdata object
save(denver.data.B, denver.model.B, est.denver.model.B, est.denver.B.cv,
     file = here::here("Results", "Denver_ST_Model_B.rdata"))

3.2.4 Prediction using the CV model

Making predictions using the CV model. Printing out the CV summary statistics as well

pred.B.cv <- predictCV(denver.model.B, est.denver.B.cv, LTA = T)
pred.B.cv.log <- predictCV(denver.model.B, est.denver.B.cv,
                           LTA = T, transform="unbiased")

head(pred.B.cv$pred.all$EX)
##            central       d_1      d_10      d_11      d_12      d_13      d_14
## 2008-12-29      NA 0.2221147 0.3406085 0.4326163 0.4042291 0.3058418 0.3641883
## 2009-01-05      NA 0.2362779 0.3500038 0.4331282 0.4063972 0.3093446 0.3606838
## 2009-01-12      NA 0.2901571 0.3954783 0.4669599 0.4364097 0.3407006 0.3848774
## 2009-01-19      NA 0.2460560 0.3517986 0.4190743 0.3967714 0.3024217 0.3391215
## 2009-01-26      NA 0.2449216 0.3470935 0.4075001 0.3868296 0.2938626 0.3226113
## 2009-02-02      NA 0.2523968 0.3500981 0.4036145 0.3821287 0.2905731 0.3107640
##                 d_15      d_16      d_17      d_18      d_19       d_2
## 2008-12-29 0.1678977 0.4593548 0.2136028 0.3918066 0.2466937 0.3625573
## 2009-01-05 0.1907173 0.4715380 0.2343124 0.3924895 0.2662171 0.3934257
## 2009-01-12 0.2479878 0.5120979 0.2890843 0.4273435 0.3281190 0.4650666
## 2009-01-19 0.2179013 0.4841166 0.2578214 0.3777103 0.2854634 0.4396387
## 2009-01-26 0.2239366 0.4875188 0.2624378 0.3649943 0.2886366 0.4599268
## 2009-02-02 0.2353514 0.4983295 0.2726721 0.3595736 0.3009704 0.4920244
##                 d_20      d_21      d_22      d_23      d_25      d_26
## 2008-12-29 0.4446601 0.2684749 0.2532801 0.4174355 0.3996886 0.3506621
## 2009-01-05 0.4439298 0.2795857 0.2700896 0.4237603 0.3972552 0.3569233
## 2009-01-12 0.4764637 0.3266903 0.3268774 0.4580907 0.4298767 0.4033645
## 2009-01-19 0.4271647 0.2844558 0.2849396 0.4231039 0.3764770 0.3515307
## 2009-01-26 0.4140016 0.2809114 0.2861979 0.4183198 0.3612660 0.3435924
## 2009-02-02 0.4083017 0.2847089 0.2960272 0.4194260 0.3541144 0.3451327
##                 d_28      d_29       d_3      d_31      d_32      d_33
## 2008-12-29 0.2409698 0.1995817 0.3692714 0.5036412 0.1741139 0.2036382
## 2009-01-05 0.2485556 0.2107143 0.3744607 0.5023219 0.1846407 0.2163498
## 2009-01-12 0.2912878 0.2575282 0.4074411 0.5412665 0.2291167 0.2716511
## 2009-01-19 0.2481723 0.2143395 0.3706515 0.4834266 0.1886138 0.2228577
## 2009-01-26 0.2435013 0.2088262 0.3633731 0.4703971 0.1842667 0.2206076
## 2009-02-02 0.2472127 0.2093470 0.3610948 0.4683630 0.1862272 0.2284429
##                 d_34      d_35      d_36      d_37      d_38      d_39
## 2008-12-29 0.4128579 0.3852355 0.3778185 0.3209316 0.3397267 0.4452858
## 2009-01-05 0.4182664 0.3859322 0.3777464 0.3242407 0.3488059 0.4394305
## 2009-01-12 0.4663383 0.4208596 0.4119223 0.3607388 0.3920069 0.4736202
## 2009-01-19 0.4104876 0.3701521 0.3617756 0.3152357 0.3505981 0.4105697
## 2009-01-26 0.4014505 0.3559547 0.3487988 0.3056072 0.3459112 0.3916376
## 2009-02-02 0.4028516 0.3482007 0.3434432 0.3031027 0.3482623 0.3828116
##                  d_4      d_40      d_41      d_42      d_43      d_44
## 2008-12-29 0.2337599 0.2196567 0.4394402 0.2474255 0.3687938 0.2565000
## 2009-01-05 0.2593379 0.2403631 0.4318852 0.2683613 0.3694133 0.2686698
## 2009-01-12 0.3592792 0.3351975 0.4576142 0.3636559 0.4042105 0.3208187
## 2009-01-19 0.2814001 0.2517219 0.4015666 0.2811304 0.3545340 0.2742532
## 2009-01-26 0.2866305 0.2506123 0.3817479 0.2817196 0.3417959 0.2709059
## 2009-02-02 0.3128546 0.2695410 0.3695163 0.3033119 0.3363801 0.2761605
##                 d_45      d_46      d_47      d_48      d_49       d_5
## 2008-12-29 0.4599630 0.2052082 0.2593854 0.4087507 0.3995698 0.3520960
## 2009-01-05 0.4534119 0.2186416 0.2746555 0.4097246 0.4002574 0.3606769
## 2009-01-12 0.4800935 0.2671311 0.3324100 0.4456410 0.4288115 0.4090039
## 2009-01-19 0.4248871 0.2295691 0.2858408 0.3952905 0.3877649 0.3594264
## 2009-01-26 0.4057384 0.2301338 0.2854623 0.3827498 0.3764931 0.3529371
## 2009-02-02 0.3939555 0.2386682 0.2947132 0.3777779 0.3705637 0.3552173
##                 d_50      d_51      d_52      d_53      d_54      d_55
## 2008-12-29 0.1861633 0.3432858 0.4579974 0.2042064 0.2805190 0.2831990
## 2009-01-05 0.1976868 0.3443546 0.4551332 0.2352984 0.2928957 0.2966881
## 2009-01-12 0.2443557 0.3795868 0.4855001 0.3410345 0.3477630 0.3843874
## 2009-01-19 0.2051665 0.3290406 0.4339684 0.2695353 0.2983259 0.2939665
## 2009-01-26 0.2044018 0.3147844 0.4184742 0.2820350 0.2951119 0.2862090
## 2009-02-02 0.2119916 0.3066881 0.4103175 0.3166709 0.3015707 0.2988774
##                 d_56      d_57      d_58       d_6      d_60      d_61
## 2008-12-29 0.3236820 0.4069372 0.2873552 0.3110044 0.3325707 0.3335435
## 2009-01-05 0.3325734 0.4028505 0.2997368 0.3172467 0.3407407 0.3441049
## 2009-01-12 0.3812601 0.4347464 0.3480880 0.3636228 0.3914413 0.3906665
## 2009-01-19 0.3321456 0.3773114 0.3070465 0.3116255 0.3379317 0.3479038
## 2009-01-26 0.3262773 0.3585939 0.3046115 0.3033735 0.3307930 0.3438547
## 2009-02-02 0.3293816 0.3472629 0.3094101 0.3044064 0.3335203 0.3471777
##                 d_62      d_63      d_64      d_65      d_66      d_67
## 2008-12-29 0.2760306 0.3336212 0.4020540 0.3403297 0.2661197 0.2669488
## 2009-01-05 0.2896380 0.3397206 0.3955047 0.3478145 0.2773346 0.2899485
## 2009-01-12 0.3431996 0.3801296 0.4220059 0.3893675 0.3283769 0.3875274
## 2009-01-19 0.2979911 0.3350599 0.3662329 0.3462004 0.2816821 0.3077500
## 2009-01-26 0.2959139 0.3267336 0.3459267 0.3395877 0.2783287 0.3117944
## 2009-02-02 0.3023260 0.3251484 0.3322248 0.3397984 0.2840694 0.3377518
##                 d_68       d_7       d_8 d_24 d_27 d_30 d_59 d_9
## 2008-12-29 0.2688478 0.3396584 0.1683511   NA   NA   NA   NA  NA
## 2009-01-05 0.2823005 0.3492668 0.1880529   NA   NA   NA   NA  NA
## 2009-01-12 0.3301128 0.3949132 0.2421639   NA   NA   NA   NA  NA
## 2009-01-19 0.2925372 0.3513175 0.2109790   NA   NA   NA   NA  NA
## 2009-01-26 0.2918368 0.3465619 0.2168041   NA   NA   NA   NA  NA
## 2009-02-02 0.2980425 0.3493424 0.2297050   NA   NA   NA   NA  NA
tail(pred.B.cv$pred.all$EX)
##            central       d_1      d_10      d_11      d_12      d_13      d_14
## 2020-11-02      NA 1.0150375 1.0776936 1.0198085 1.0591903 0.9841650 0.9620116
## 2020-11-09      NA 0.6251803 0.6913531 0.6313462 0.6727281 0.5863371 0.5937758
## 2020-11-16      NA 1.0243758 1.1019127 1.0535331 1.0886497 1.0044736 1.0171267
## 2020-11-23      NA 0.9527293 1.0288895 0.9916107 1.0192567 0.9329078 0.9577483
## 2020-11-30      NA 1.1483583 1.2294370 1.2083545 1.2296129 1.1446814 1.1793494
## 2020-12-07      NA        NA        NA        NA        NA        NA        NA
##                 d_15     d_16     d_17     d_18     d_19       d_2      d_20
## 2020-11-02 1.1561735 1.147123 1.164793 1.020855 1.116269 1.2825994 1.0213387
## 2020-11-09 0.7553162 0.723739 0.758770 0.643016 0.725545 0.8196408 0.6443358
## 2020-11-16 1.1506762 1.125905 1.152468 1.067550 1.125433 1.1909961 1.0745858
## 2020-11-23 1.0582438 1.017714 1.065531 1.007662 1.043340 1.0731742 1.0132734
## 2020-11-30 1.2396615 1.214072 1.253733 1.224595 1.228521 1.2290675 1.2376935
## 2020-12-07        NA       NA       NA       NA       NA        NA        NA
##                 d_21      d_22      d_23      d_25      d_26      d_28
## 2020-11-02 1.0448420 1.0920563 1.0911951 0.9647147 1.0342517 0.9003043
## 2020-11-09 0.6539470 0.6925186 0.6945855 0.5946751 0.6377978 0.5131376
## 2020-11-16 1.0615172 1.0932745 1.1080573 1.0242741 1.0620301 0.9198463
## 2020-11-23 0.9891587 1.0159927 1.0238217 0.9666114 1.0002494 0.8438110
## 2020-11-30 1.1895191 1.2100321 1.2224294 1.1876949 1.2141364 1.0445219
## 2020-12-07        NA        NA        NA        NA        NA        NA
##                 d_29       d_3      d_31      d_32      d_33      d_34
## 2020-11-02 1.0361800 1.0783991 1.0059734 0.9684093 0.9456605 1.0617426
## 2020-11-09 0.6552218 0.6833334 0.6340943 0.5773669 0.5591726 0.6839008
## 2020-11-16 1.0648836 1.0983172 1.0610828 0.9879177 0.9621700 1.1041349
## 2020-11-23 0.9970938 1.0236954 1.0015543 0.9135634 0.8921258 1.0434900
## 2020-11-30 1.2023364 1.2309221 1.2138238 1.1171841 1.0893082 1.2555070
## 2020-12-07        NA        NA        NA        NA        NA        NA
##                 d_35      d_36      d_37      d_38      d_39       d_4
## 2020-11-02 1.0620075 0.9817250 0.9780876 1.0803897 0.9236184 1.0462908
## 2020-11-09 0.6892449 0.6046740 0.5928632 0.6885219 0.5715089 0.6454975
## 2020-11-16 1.1178506 1.0294449 1.0136770 1.0997828 0.9983757 1.0551971
## 2020-11-23 1.0615771 0.9690825 0.9455575 1.0237380 0.9552188 0.9971032
## 2020-11-30 1.2839672 1.1857704 1.1568365 1.2283501 1.1776324 1.1925871
## 2020-12-07        NA        NA        NA        NA        NA        NA
##                 d_40      d_41      d_42      d_43      d_44      d_45
## 2020-11-02 0.9054965 0.9479669 1.0124860 0.9970620 1.0414596 0.9870247
## 2020-11-09 0.5470347 0.5803466 0.6255479 0.6206986 0.6487289 0.6205416
## 2020-11-16 0.9416465 1.0204236 1.0446210 1.0461720 1.0589809 1.0570829
## 2020-11-23 0.8801214 0.9676142 0.9917533 0.9838337 0.9898298 1.0052004
## 2020-11-30 1.0765392 1.1988128 1.1928202 1.2017317 1.1926276 1.2333951
## 2020-12-07        NA        NA        NA        NA        NA        NA
##                 d_46      d_47      d_48      d_49       d_5      d_50
## 2020-11-02 0.9517694 1.0548802 1.0339263 1.0387043 1.0456600 0.8876971
## 2020-11-09 0.5529526 0.6716038 0.6595257 0.6528972 0.6698027 0.4947756
## 2020-11-16 0.9489201 1.0755656 1.0816021 1.0767264 1.0721795 0.8906770
## 2020-11-23 0.8679939 1.0014233 1.0234128 1.0070041 1.0084142 0.8111266
## 2020-11-30 1.0595220 1.1872085 1.2388418 1.2182878 1.2118826 1.0028663
## 2020-12-07        NA        NA        NA        NA        NA        NA
##                 d_51      d_52      d_53      d_54      d_55      d_56
## 2020-11-02 1.0327922 1.0314674 1.1132033 1.0377687 0.9845357 1.0231886
## 2020-11-09 0.6669587 0.6537372 0.6812753 0.6531501 0.6206208 0.6420580
## 2020-11-16 1.0959861 1.0842459 1.0826501 1.0651263 1.0527421 1.0461825
## 2020-11-23 1.0403651 1.0286754 1.0082953 0.9966882 1.0133027 0.9831596
## 2020-11-30 1.2553539 1.2509919 1.1861315 1.2032899 1.2285154 1.1868197
## 2020-12-07        NA        NA        NA        NA        NA        NA
##                 d_57      d_58       d_6      d_60      d_61      d_62
## 2020-11-02 0.9935330 1.0948235 0.9957890 1.0276878 1.0841772 1.0625656
## 2020-11-09 0.6373275 0.6998669 0.6120707 0.6484534 0.6950927 0.6726742
## 2020-11-16 1.0673926 1.1071086 1.0330153 1.0615123 1.0978030 1.0729691
## 2020-11-23 1.0219978 1.0326973 0.9690720 0.9982060 1.0288889 1.0037612
## 2020-11-30 1.2476433 1.2302258 1.1791473 1.2061263 1.2307727 1.2029491
## 2020-12-07        NA        NA        NA        NA        NA        NA
##                 d_63      d_64      d_65      d_66      d_67      d_68
## 2020-11-02 1.0552428 0.9614875 1.0744089 1.0073493 1.0145668 1.0898587
## 2020-11-09 0.6727983 0.6009151 0.6723145 0.6174771 0.6280493 0.6869631
## 2020-11-16 1.0878553 1.0364567 1.0906425 1.0263662 1.0443520 1.0910055
## 2020-11-23 1.0226460 0.9897999 1.0246298 0.9530366 0.9796384 1.0119553
## 2020-11-30 1.2316994 1.2194697 1.2377665 1.1550219 1.1658121 1.2070013
## 2020-12-07        NA        NA        NA        NA        NA        NA
##                  d_7       d_8 d_24 d_27 d_30 d_59 d_9
## 2020-11-02 1.0722112 1.0336011   NA   NA   NA   NA  NA
## 2020-11-09 0.6829798 0.6237523   NA   NA   NA   NA  NA
## 2020-11-16 1.0884182 1.0135626   NA   NA   NA   NA  NA
## 2020-11-23 1.0188960 0.9224215   NA   NA   NA   NA  NA
## 2020-11-30 1.2223381 1.1060283   NA   NA   NA   NA  NA
## 2020-12-07        NA        NA   NA   NA   NA   NA  NA
head(pred.B.cv.log$pred.all$EX)
##            central      d_1     d_10     d_11     d_12     d_13     d_14
## 2008-12-29      NA 1.282235 1.442799 1.583269 1.539445 1.395196 1.479022
## 2009-01-05      NA 1.300224 1.456188 1.583662 1.542359 1.399704 1.473440
## 2009-01-12      NA 1.371921 1.523724 1.637776 1.588965 1.443938 1.509156
## 2009-01-19      NA 1.312506 1.458428 1.560892 1.526900 1.389425 1.441364
## 2009-01-26      NA 1.310834 1.451441 1.542684 1.511546 1.377356 1.417528
## 2009-02-02      NA 1.320527 1.455698 1.536513 1.504268 1.372661 1.400658
##                d_15     d_16     d_17     d_18     d_19      d_2     d_20
## 2008-12-29 1.213964 1.626691 1.271460 1.519599 1.313545 1.476867 1.602453
## 2009-01-05 1.241845 1.646174 1.297752 1.520258 1.339146 1.522727 1.600861
## 2009-01-12 1.314912 1.713899 1.370525 1.573836 1.424383 1.635405 1.653416
## 2009-01-19 1.275833 1.666266 1.328104 1.497353 1.364677 1.594005 1.573571
## 2009-01-26 1.283467 1.671668 1.334058 1.478216 1.368828 1.626395 1.552746
## 2009-02-02 1.298132 1.689626 1.347634 1.470061 1.385670 1.679228 1.543732
##                d_21     d_22     d_23     d_25     d_26     d_28     d_29
## 2008-12-29 1.342389 1.323985 1.559910 1.530633 1.459404 1.305992 1.253021
## 2009-01-05 1.357173 1.346040 1.569373 1.526742 1.468146 1.315789 1.266849
## 2009-01-12 1.422434 1.424330 1.623792 1.577212 1.537547 1.373098 1.327381
## 2009-01-19 1.363448 1.365540 1.567641 1.495072 1.459569 1.315044 1.271122
## 2009-01-26 1.358491 1.367026 1.559902 1.472401 1.447781 1.308825 1.264010
## 2009-02-02 1.363556 1.380350 1.561433 1.461831 1.449824 1.313622 1.264573
##                 d_3     d_31     d_32     d_33     d_34     d_35     d_36
## 2008-12-29 1.486559 1.699155 1.222324 1.258190 1.550991 1.509530 1.498491
## 2009-01-05 1.493880 1.696522 1.234951 1.274004 1.559057 1.510217 1.498009
## 2009-01-12 1.543595 1.763537 1.290835 1.346180 1.635516 1.563566 1.549752
## 2009-01-19 1.487534 1.664141 1.239369 1.281859 1.546419 1.485993 1.473682
## 2009-01-26 1.476503 1.642368 1.233811 1.278805 1.532299 1.464834 1.454468
## 2009-02-02 1.472958 1.638853 1.236095 1.288729 1.534286 1.453361 1.446538
##                d_37     d_38     d_39      d_4     d_40     d_41     d_42
## 2008-12-29 1.415959 1.442484 1.602837 1.297257 1.279090 1.594110 1.315106
## 2009-01-05 1.420277 1.455278 1.593111 1.330500 1.305492 1.581694 1.342560
## 2009-01-12 1.472731 1.519193 1.648186 1.469989 1.435012 1.622541 1.476437
## 2009-01-19 1.406942 1.457302 1.547208 1.359575 1.319818 1.533802 1.359208
## 2009-01-26 1.393238 1.450274 1.517979 1.366480 1.318138 1.503463 1.359786
## 2009-02-02 1.389583 1.453526 1.504477 1.402613 1.343157 1.485004 1.389292
##                d_43     d_44     d_45     d_46     d_47     d_48     d_49
## 2008-12-29 1.485028 1.328255 1.627164 1.260113 1.330323 1.544567 1.532289
## 2009-01-05 1.485578 1.344130 1.616112 1.277011 1.350494 1.545898 1.532918
## 2009-01-12 1.537846 1.415726 1.659428 1.340326 1.430508 1.602272 1.576937
## 2009-01-19 1.463049 1.351026 1.569991 1.290806 1.365192 1.523466 1.513210
## 2009-01-26 1.444318 1.346280 1.539968 1.291446 1.364490 1.504376 1.496002
## 2009-02-02 1.436357 1.353197 1.521743 1.302445 1.377027 1.496835 1.486972
##                 d_5     d_50     d_51     d_52     d_53     d_54     d_55
## 2008-12-29 1.460218 1.236341 1.447515 1.623969 1.259479 1.358736 1.363004
## 2009-01-05 1.472461 1.250530 1.448713 1.618896 1.298897 1.375353 1.381134
## 2009-01-12 1.545054 1.310145 1.500347 1.668424 1.443413 1.452640 1.507365
## 2009-01-19 1.470068 1.259688 1.426140 1.584314 1.343539 1.382343 1.376768
## 2009-01-26 1.460354 1.258638 1.405751 1.559706 1.360215 1.377721 1.365904
## 2009-02-02 1.463530 1.268160 1.394264 1.546847 1.407976 1.386502 1.383145
##                d_56     d_57     d_58      d_6     d_60     d_61     d_62
## 2008-12-29 1.419311 1.541743 1.367975 1.402660 1.431334 1.432641 1.354452
## 2009-01-05 1.431656 1.535212 1.384799 1.411036 1.442756 1.447623 1.372612
## 2009-01-12 1.502778 1.584749 1.453199 1.477641 1.517495 1.516410 1.447769
## 2009-01-19 1.430506 1.496115 1.394599 1.402472 1.438191 1.452759 1.383480
## 2009-01-26 1.421936 1.468229 1.391072 1.390708 1.427767 1.446747 1.380372
## 2009-02-02 1.426203 1.451576 1.397657 1.391965 1.431516 1.451453 1.389072
##                d_63     d_64     d_65     d_66     d_67     d_68      d_7
## 2008-12-29 1.433593 1.535613 1.443354 1.339919 1.341034 1.343678 1.441428
## 2009-01-05 1.442016 1.525186 1.453836 1.354718 1.371857 1.361547 1.455115
## 2009-01-12 1.501162 1.565782 1.515189 1.425371 1.512106 1.427925 1.522864
## 2009-01-19 1.434750 1.480553 1.450907 1.360109 1.395876 1.375020 1.457727
## 2009-01-26 1.422649 1.450561 1.441132 1.355365 1.401302 1.373860 1.450669
## 2009-02-02 1.420241 1.430646 1.441276 1.363020 1.437972 1.382261 1.454598
##                 d_8 d_24 d_27 d_30 d_59 d_9
## 2008-12-29 1.215206   NA   NA   NA   NA  NA
## 2009-01-05 1.239086   NA   NA   NA   NA  NA
## 2009-01-12 1.307705   NA   NA   NA   NA  NA
## 2009-01-19 1.267327   NA   NA   NA   NA  NA
## 2009-01-26 1.274548   NA   NA   NA   NA  NA
## 2009-02-02 1.290957   NA   NA   NA   NA  NA
tail(pred.B.cv.log$pred.all$EX)
##            central      d_1     d_10     d_11     d_12     d_13     d_14
## 2020-11-02      NA 2.774707 2.950054 2.785553 2.898576 2.687846 2.630276
## 2020-11-09      NA 1.878626 2.004501 1.888545 1.969049 1.805260 1.819653
## 2020-11-16      NA 2.800152 3.022040 2.880341 2.984271 2.742101 2.778414
## 2020-11-23      NA 2.606599 2.809342 2.707406 2.784120 2.552640 2.618159
## 2020-11-30      NA 3.170127 3.433569 3.362966 3.436103 3.154885 3.267818
## 2020-12-07      NA       NA       NA       NA       NA       NA       NA
##                d_15     d_16     d_17     d_18     d_19      d_2     d_20
## 2020-11-02 3.190019 3.160998 3.220468 2.788032 3.065765 3.621659 2.788091
## 2020-11-09 2.136362 2.069476 2.145388 1.910367 2.073920 2.279030 1.912050
## 2020-11-16 3.172268 3.093627 3.180104 2.920361 3.093432 3.303470 2.939799
## 2020-11-23 2.892275 2.776311 2.915228 2.750464 2.849728 2.936171 2.764980
## 2020-11-30 3.467849 3.378837 3.519052 3.416862 3.429857 3.431631 3.460948
## 2020-12-07       NA       NA       NA       NA       NA       NA       NA
##                d_21     d_22     d_23     d_25     d_26     d_28     d_29
## 2020-11-02 2.855439 2.997261 2.989548 2.635053 2.826848 2.471497 2.831633
## 2020-11-09 1.931396 2.009606 2.010341 1.819930 1.901200 1.677983 1.934420
## 2020-11-16 2.903134 2.999853 3.039397 2.796531 2.905443 2.520062 2.913769
## 2020-11-23 2.700599 2.776625 2.793773 2.639914 2.731252 2.335620 2.722903
## 2020-11-30 3.300045 3.371342 3.407747 3.293342 3.382729 2.854969 3.343584
## 2020-12-07       NA       NA       NA       NA       NA       NA       NA
##                 d_3     d_31     d_32     d_33     d_34     d_35     d_36
## 2020-11-02 2.953556 2.746222 2.645493 2.587066 2.903588 2.905010 2.680959
## 2020-11-09 1.989208 1.893068 1.788923 1.757522 1.989679 2.000687 1.838449
## 2020-11-16 3.011998 2.901200 2.696736 2.629663 3.028781 3.070965 2.811085
## 2020-11-23 2.795333 2.733588 2.503372 2.451859 2.850661 2.902839 2.646291
## 2020-11-30 3.439163 3.380359 3.068778 2.986619 3.524298 3.625967 3.286644
## 2020-12-07       NA       NA       NA       NA       NA       NA       NA
##                d_37     d_38     d_39      d_4     d_40     d_41     d_42
## 2020-11-02 2.670951 2.958414 2.532819 2.862185 2.485235 2.591327 2.766262
## 2020-11-09 1.816719 1.998874 1.780825 1.916678 1.736214 1.793864 1.878288
## 2020-11-16 2.766993 3.015369 2.728852 2.886909 2.575932 2.785324 2.855727
## 2020-11-23 2.584795 2.794431 2.613638 2.723922 2.422177 2.642061 2.708619
## 2020-11-30 3.193170 3.428972 3.264977 3.312219 2.948057 3.329583 3.312052
## 2020-12-07       NA       NA       NA       NA       NA       NA       NA
##                d_43     d_44     d_45     d_46     d_47     d_48     d_49
## 2020-11-02 2.722106 2.849035 2.695286 2.602890 2.883439 2.825016 2.837492
## 2020-11-09 1.867950 1.923271 1.867953 1.746722 1.965163 1.942641 1.928815
## 2020-11-16 2.858200 2.898368 2.890126 2.595271 2.943178 2.962719 2.946500
## 2020-11-23 2.685333 2.704593 2.744014 2.393588 2.732954 2.795321 2.747977
## 2020-11-30 3.339171 3.312770 3.447694 2.899082 3.291298 3.467551 3.394645
## 2020-12-07       NA       NA       NA       NA       NA       NA       NA
##                 d_5     d_50     d_51     d_52     d_53     d_54     d_55
## 2020-11-02 2.861379 2.441746 2.820986 2.817935 3.056913 2.834571 2.690471
## 2020-11-09 1.964621 1.648271 1.956328 1.931113 1.984325 1.929266 1.869373
## 2020-11-16 2.937661 2.448831 3.004143 2.969877 2.964019 2.912666 2.879503
## 2020-11-23 2.756239 2.261637 2.841529 2.809354 2.751570 2.720093 2.768092
## 2020-11-30 3.378500 2.739844 3.523211 3.509101 3.287312 3.344715 3.432994
## 2020-12-07       NA       NA       NA       NA       NA       NA       NA
##                d_56     d_57     d_58      d_6     d_60     d_61     d_62
## 2020-11-02 2.798000 2.714004 3.000952 2.720691 2.806878 2.971894 2.911998
## 2020-11-09 1.911002 1.900529 2.021591 1.853257 1.920733 2.013808 1.971365
## 2020-11-16 2.862483 2.921727 3.037711 2.822880 2.902924 3.012335 2.941412
## 2020-11-23 2.687699 2.792176 2.819993 2.647907 2.724943 2.811850 2.744603
## 2020-11-30 3.295118 3.499320 3.436196 3.267027 3.355099 3.441227 3.349662
## 2020-12-07       NA       NA       NA       NA       NA       NA       NA
##                d_63     d_64     d_65     d_66     d_67     d_68      d_7
## 2020-11-02 2.886553 2.628531 2.941277 2.750721 2.769566 2.986962 2.935527
## 2020-11-09 1.968822 1.832489 1.967076 1.862357 1.881321 1.996064 1.988873
## 2020-11-16 2.981386 2.832422 2.988446 2.802941 2.852425 2.989528 2.983163
## 2020-11-23 2.793093 2.703318 2.797404 2.604808 2.673627 2.762224 2.782926
## 2020-11-30 3.442660 3.401574 3.462007 3.188151 3.220937 3.357255 3.411142
## 2020-12-07       NA       NA       NA       NA       NA       NA       NA
##                 d_8 d_24 d_27 d_30 d_59 d_9
## 2020-11-02 2.824528   NA   NA   NA   NA  NA
## 2020-11-09 1.874438   NA   NA   NA   NA  NA
## 2020-11-16 2.767696   NA   NA   NA   NA  NA
## 2020-11-23 2.526525   NA   NA   NA   NA  NA
## 2020-11-30 3.035855   NA   NA   NA   NA  NA
## 2020-12-07       NA   NA   NA   NA   NA  NA
summary(pred.B.cv)
## Cross-validation predictions for STmodel with 10 CV-groups.
##   Predictions for 790 observations.
##   Temporal averages for 63 locations.
## 
## RMSE:
##              EX.mu EX.mu.beta         EX
## obs     0.13784779 0.13784779 0.11072223
## average 0.08626559 0.08626559 0.06637898
## 
## R2:
##             EX.mu EX.mu.beta        EX
## obs     0.7179091  0.7179091 0.8180051
## average 0.6694103  0.6694103 0.8042619
## 
## Coverage of 95% prediction intervals:
##                EX
## obs     0.9367089
## average 0.8571429
summary(pred.B.cv.log)
## Cross-validation predictions for STmodel with 10 CV-groups.
##   Predictions for 790 observations.
##   Temporal averages for 63 locations.
## 
## RMSE:
##             EX.mu EX.mu.beta         EX    EX.pred
## obs     0.1938374  0.1938374 0.14992738 0.14995756
## average 0.1017864  0.1017864 0.06827081 0.06883603
## 
## R2:
##             EX.mu EX.mu.beta        EX   EX.pred
## obs     0.6931546  0.6931546 0.8164280 0.8163541
## average 0.6921204  0.6921204 0.8614931 0.8591901
## 
## Coverage of 95% prediction intervals:
##           EX.pred
## obs     0.9367089
## average 0.8730159
par(mfrow=c(1,2), mar=c(3.3,3.3,1.5,1), mgp=c(2,1,0))
plot(pred.B.cv, "obs", ID="all", pch=c(19,NA), cex=.25, lty=c(NA,2),
     col=c("ID", "black", "grey"),
     ylim=c(-1,2),
     xlab="Observations", ylab="Predictions",
     main="Cross-validation BC (log ug/m3)")
with(pred.B.cv.log$pred.LTA, plotCI(obs, EX.pred, uiw=1.96*sqrt(VX.pred),
                                    xlab="Observations", ylab="Predictions",
                                    main="Temporal average BC (ug/m3)"))
abline(0, 1, col="grey")

jpeg(filename = here::here("Figs", "ST_CV_Obs_vs_Pred_BC_ModB.jpeg"),
     width = 8, height = 4, units = "in", res = 500)
par(mfrow=c(1,2), mar=c(3.3,3.3,1.5,1), mgp=c(2,1,0))
plot(pred.B.cv, "obs", ID="all", pch=c(19,NA), cex=.25, lty=c(NA,2),
     col=c("ID", "black", "grey"),
     ylim=c(-1,2),
     xlab="Observations", ylab="Predictions",
     main="Cross-validation BC (log ug/m3)")
with(pred.B.cv.log$pred.LTA, plotCI(obs, EX.pred, uiw=1.96*sqrt(VX.pred),
                                      xlab="Observations", ylab="Predictions",
                                      main="Temporal average BC (ug/m3)"))
abline(0, 1, col="grey")
dev.off()
## quartz_off_screen 
##                 2

3.2.5 Predictions for 2009-2020

What do the long-term predictions look like for this model? Predicting at the distributed (residential + community) sites

x.B <- coef(est.denver.model.B, pars = "cov")$par
x.B
## [1] -10.463290  -7.028546  -7.730966  13.512563  -2.978499  -4.934033
E.B <- predict(denver.model.B, est.denver.model.B, LTA = T,
               transform="unbiased", pred.var=FALSE)
print(E.B)
## Prediction for STmodel.
## 
## Regression parameters:
##   0 Spatio-temporal covariate(s).
##   42 beta-fields regression parameters in x$pars.
## 
## Regression parameters are assumed to be known and
##  prediction variances do NOT include
##  uncertainties in regression parameters.
## 
## Prediction of beta-fields, (x$beta):
## List of 2
##  $ mu: num [1:69, 1:3] 0.403 0.203 0.276 0.237 0.252 ...
##   ..- attr(*, "dimnames")=List of 2
##  $ EX: num [1:69, 1:3] 0.404 0.204 0.275 0.236 0.25 ...
##   ..- attr(*, "dimnames")=List of 2
## 
## Predictions for log-Gaussian field of type: unbiased 
## 
## Predictions for 622 times at 69 locations.
## List of 5
##  $ EX.mu     : num [1:622, 1:69] 1.5 1.53 1.62 1.56 1.58 ...
##   ..- attr(*, "dimnames")=List of 2
##  $ EX.mu.beta: num [1:622, 1:69] 1.54 1.57 1.65 1.6 1.61 ...
##   ..- attr(*, "dimnames")=List of 2
##  $ EX        : num [1:622, 1:69] 1.58 1.61 1.7 1.64 1.65 ...
##   ..- attr(*, "dimnames")=List of 2
##  $ EX.pred   : num [1:622, 1:69] 1.59 1.61 1.7 1.64 1.65 ...
##   ..- attr(*, "dimnames")=List of 2
##  $ log.EX    : num [1:622, 1:69] 0.432 0.448 0.503 0.467 0.474 ...
##   ..- attr(*, "dimnames")=List of 2
## 
## Variances have been computed.
## List of 2
##  $ log.VX     : num [1:622, 1:69] 0.0514 0.0513 0.0512 0.0511 0.0511 ...
##   ..- attr(*, "dimnames")=List of 2
##  $ log.VX.pred: num [1:622, 1:69] 0.0586 0.0585 0.0584 0.0583 0.0583 ...
##   ..- attr(*, "dimnames")=List of 2
## 
## Mean squared prediciton errors NOT been computed.
## 
## 69 temporal averages have been compute.
## List of 1
##  $ LTA:'data.frame': 69 obs. of  4 variables:
week_preds <- as.data.frame(E.B$EX) %>%
  mutate(week = as.Date(rownames(E.B$EX)),
         year = year(as.Date(rownames(E.B$EX))))

week_sites <- colnames(week_preds)[which(str_detect(colnames(week_preds), "d_"))]

box_preds <- week_preds %>%
  select(week, all_of(week_sites)) %>%
  #filter(week %in% week_weeks) %>%
  mutate(month = month(week),
         year = year(week)) %>%
  filter(year %in% c(2009:2020))

box_data <- box_preds %>%
  pivot_longer(-c(week, month, year), names_to = "location", values_to = "pred")
summary(box_data)
##       week                month             year        location        
##  Min.   :2009-01-05   Min.   : 1.000   Min.   :2009   Length:42228      
##  1st Qu.:2012-01-09   1st Qu.: 4.000   1st Qu.:2012   Class :character  
##  Median :2014-12-29   Median : 7.000   Median :2014   Mode  :character  
##  Mean   :2014-12-28   Mean   : 6.494   Mean   :2014                     
##  3rd Qu.:2017-12-18   3rd Qu.: 9.000   3rd Qu.:2017                     
##  Max.   :2020-12-07   Max.   :12.000   Max.   :2020                     
##                                                                         
##       pred       
##  Min.   :0.3553  
##  1st Qu.:1.1765  
##  Median :1.4744  
##  Mean   :1.5943  
##  3rd Qu.:2.0065  
##  Max.   :4.0346  
##  NA's   :861
#' Box plot of weekly BC predicted at all distributed sites grouped by month
box_summary <- box_data %>%
  group_by(month) %>%
  summarize(median = median(pred, na.rm=T)) %>%
  arrange(desc(median))
box_summary
pred_box_plot <- ggplot(box_data) +
  geom_boxplot(aes(x = as.factor(month), y = pred), #, color = as.factor(month)),
               show.legend = F) +
  #scale_color_viridis(name = "Month", discrete = T) +
  facet_wrap(. ~ year, scales = "free") +
  xlab("") + ylab("Distributed site BC concentration (\u03BCg/m\u00B3): Predicted") +
  # theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust = 1)) +
  scale_x_discrete(labels = str_sub(month.name, 1, 1))
pred_box_plot

4 Summary Table

Summary of model diagnostics
model trend_data cov_beta0 cov_beta1 cov_beta2 cov_nu_error no_basis_fns df_per_year all_converge cv_rmse_obs cv_rmse_avg cv_r2_obs cv_r2_avg cv_coverage_obs cv_coverage_avg
Model A BC + NO2 + PM2.5 iid iid NA exp 1 4 TRUE 0.15 0.06 0.81 0.89 0.94 0.90
Model B BC + NO2 + PM2.5 iid iid iid exp 2 4 TRUE 0.15 0.07 0.82 0.86 0.94 0.87